SCIP-SDP  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
sdpisolver_dsdp.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of SCIPSDP - a solving framework for mixed-integer */
4 /* semidefinite programms based on SCIP. */
5 /* */
6 /* Copyright (C) 2011-2013 Discrete Optimization, TU Darmstadt */
7 /* EDOM, FAU Erlangen-Nürnberg */
8 /* 2014-2015 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-2015 Zuse Institute Berlin */
28 /* SCIP is distributed under the terms of the SCIP Academic Licence, */
29 /* see file COPYING in the SCIP distribution. */
30 /* */
31 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
32 
33 /*#define SCIP_DEBUG*/
34 /*#define SCIP_MORE_DEBUG*/
35 
41 #include <assert.h>
42 
43 #include "sdpi/sdpisolver.h"
44 
45 /* turn off warning for DSDSP */
46 #pragma GCC diagnostic ignored "-Wstrict-prototypes"
47 #include "dsdp5.h" /* for DSDPUsePenalty, etc */
48 #pragma GCC diagnostic warning "-Wstrict-prototypes"
49 
50 
51 #include "blockmemshell/memory.h" /* for memory allocation */
52 #include "scip/def.h" /* for SCIP_Real, _Bool, ... */
53 #include "scip/pub_misc.h" /* for sorting */
54 
56 #define DSDP_CALL(x) do \
57  { \
58  int _dsdperrorcode_; \
59  if ( (_dsdperrorcode_ = (x)) != 0 ) \
60  { \
61  SCIPerrorMessage("DSDP-Error <%d> in function call.\n", _dsdperrorcode_); \
62  return SCIP_LPERROR; \
63  } \
64  } \
65  while( FALSE )
66 
68 #define DSDP_CALL_BOOL(x) do \
69  { \
70  int _dsdperrorcode_; \
71  if ( (_dsdperrorcode_ = (x)) != 0 ) \
72  { \
73  SCIPerrorMessage("DSDP-Error <%d> in function call.\n", _dsdperrorcode_); \
74  return FALSE; \
75  } \
76  } \
77  while( FALSE )
78 
80 #define DSDP_CALLM(x) do \
81  { \
82  int _dsdperrorcode_; \
83  if ( (_dsdperrorcode_ = (x)) != 0 ) \
84  { \
85  SCIPerrorMessage("DSDP-Error <%d> in function call.\n", _dsdperrorcode_); \
86  return SCIP_NOMEMORY; \
87  } \
88  } \
89  while( FALSE )
90 
92 #define BMS_CALL(x) do \
93  { \
94  if( NULL == (x) ) \
95  { \
96  SCIPerrorMessage("No memory in function call.\n"); \
97  return SCIP_NOMEMORY; \
98  } \
99  } \
100  while( FALSE )
101 
103 #define CHECK_IF_SOLVED(sdpisolver) do \
104  { \
105  if (!(sdpisolver->solved)) \
106  { \
107  SCIPerrorMessage("Tried to access solution information for SDP %d ahead of solving!\n", sdpisolver->sdpcounter); \
108  return SCIP_LPERROR; \
109  } \
110  } \
111  while( FALSE )
112 
114 #define CHECK_IF_SOLVED_BOOL(sdpisolver) do \
115  { \
116  if (!(sdpisolver->solved)) \
117  { \
118  SCIPerrorMessage("Tried to access solution information for SDP %d ahead of solving!\n", sdpisolver->sdpcounter); \
119  return FALSE; \
120  } \
121  } \
122  while( FALSE )
123 
124 
126 struct SCIP_SDPiSolver
127 {
128  SCIP_MESSAGEHDLR* messagehdlr;
129  BMS_BLKMEM* blkmem;
130  DSDP dsdp;
131  SDPCone sdpcone;
132  LPCone lpcone;
133  BCone bcone;
134  int nvars;
135  int nactivevars;
136  int* inputtodsdpmapper;
139  int* dsdptoinputmapper;
140  SCIP_Real* fixedvarsval;
141  SCIP_Real fixedvarsobjcontr;
142  SCIP_Bool solved;
143  int sdpcounter;
144  SCIP_Real epsilon;
145  SCIP_Real feastol;
146  SCIP_Real objlimit;
147  SCIP_Bool sdpinfo;
148  SCIP_Bool penaltyworbound;
149 };
150 
151 
152 /*
153  * Local Functions
154  */
155 
159 static
161  int i,
162  int j
163  )
164 {
165  assert( j >= 0 );
166  assert( i >= j );
167 
168  return i*(i+1)/2 + j;
169 }
170 
171 #ifndef NDEBUG
172 
173 static
174 SCIP_Bool isFixed(
175  SCIP_SDPISOLVER* sdpisolver,
176  SCIP_Real lb,
177  SCIP_Real ub
178  )
179 {
180  assert( lb < ub + sdpisolver->feastol );
181 
182  return (REALABS(ub-lb) <= sdpisolver->feastol);
183 }
184 #else
185 #define isFixed(sdpisolver,lb,ub) (REALABS(ub-lb) <= sdpisolver->feastol)
186 #endif
187 
189 static
191  int* row, /* row indices */
192  int* col, /* column indices */
193  SCIP_Real* val, /* values */
194  int length /* length of the given arrays */
195  )
196 {
197  int firstentry;
198  int nextentry = 0;
199 
200  /* first sort by col indices */
201  SCIPsortIntIntReal(col, row, val, length);
202 
203  /* for those with identical col-indices now sort by non-decreasing row-index, first find all entries with the same col-index */
204  while (nextentry < length)
205  {
206  firstentry = nextentry; /* the next col starts where the last one ended */
207 
208  while (nextentry < length && col[nextentry] == col[firstentry]) /* as long as the row still matches, increase nextentry */
209  ++nextentry;
210 
211  /* now sort all entries between firstentry and nextentry-1 by their row-indices */
212  SCIPsortIntReal(row + firstentry, val + firstentry, nextentry - firstentry);
213  }
214 }
215 
216 
217 /*
218  * Miscellaneous Methods
219  */
220 
227  void
228  )
229 {
230  return "DSDP"; /* getting the version is not supported in DSDP */
231 }
232 
235  void
236  )
237 {
238  return "Dual-Scaling Interior Point Solver for Semidefinite Programming developed by Steve Benson, Yinyu Ye, and Xiong Zhang "
239  "(http://www.mcs.anl.gov/hs/software/DSDP/)";
240 }
241 
249  SCIP_SDPISOLVER* sdpisolver
250  )
251 {
252  assert( sdpisolver != NULL );
253  return (void*) sdpisolver->dsdp;
254 }
255 
259 /*
260  * SDPI Creation and Destruction Methods
261  */
262 
267 SCIP_RETCODE SCIPsdpiSolverCreate(
268  SCIP_SDPISOLVER** sdpisolver,
269  SCIP_MESSAGEHDLR* messagehdlr,
270  BMS_BLKMEM* blkmem
271  )
272 {
273  assert( sdpisolver != NULL );
274  assert( blkmem != NULL );
275 
276  SCIPdebugMessage("Calling SCIPsdpiCreate \n");
277 
278  BMS_CALL( BMSallocBlockMemory(blkmem, sdpisolver) );
279 
280  (*sdpisolver)->messagehdlr = messagehdlr;
281  (*sdpisolver)->blkmem = blkmem;
282 
283  /* the following four variables will be properly initialized only immediatly prior to solving because DSDP and the
284  * SDPCone need information about the number of variables and sdpblocks during creation */
285  (*sdpisolver)->dsdp = NULL;
286  (*sdpisolver)->sdpcone = NULL;
287  (*sdpisolver)->lpcone = NULL;
288  (*sdpisolver)->bcone = NULL;
289 
290  (*sdpisolver)->nvars = 0;
291  (*sdpisolver)->nactivevars = 0;
292  (*sdpisolver)->inputtodsdpmapper = NULL;
293  (*sdpisolver)->dsdptoinputmapper = NULL;
294  (*sdpisolver)->fixedvarsval = NULL;
295  (*sdpisolver)->fixedvarsobjcontr = 0.0;
296  (*sdpisolver)->solved = FALSE;
297  (*sdpisolver)->sdpcounter = 0;
298 
299  (*sdpisolver)->epsilon = 1e-5;
300  (*sdpisolver)->feastol = 1e-4;
301  (*sdpisolver)->objlimit = SCIPsdpiSolverInfinity(*sdpisolver);
302  (*sdpisolver)->sdpinfo = FALSE;
303 
304  return SCIP_OKAY;
305 }
306 
308 SCIP_RETCODE SCIPsdpiSolverFree(
309  SCIP_SDPISOLVER** sdpisolver
310  )
311 {
312  assert( sdpisolver != NULL );
313  assert( *sdpisolver != NULL );
314 
315  SCIPdebugMessage("Freeing SDPISolver\n");
316 
317  if ( (*sdpisolver)->dsdp != NULL )
318  {
319  DSDP_CALL( DSDPDestroy((*sdpisolver)->dsdp) );
320  }
321 
322  if ( (*sdpisolver)->nvars > 0 )
323  BMSfreeBlockMemoryArray((*sdpisolver)->blkmem, &(*sdpisolver)->inputtodsdpmapper, (*sdpisolver)->nvars);
324 
325  if ( (*sdpisolver)->nactivevars > 0 )
326  BMSfreeBlockMemoryArray((*sdpisolver)->blkmem, &(*sdpisolver)->dsdptoinputmapper, (*sdpisolver)->nactivevars);
327 
328  if ( (*sdpisolver)->nvars >= (*sdpisolver)->nactivevars )
329  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->fixedvarsval, (*sdpisolver)->nvars - (*sdpisolver)->nactivevars); /*lint !e776*/
330 
331  BMSfreeBlockMemory((*sdpisolver)->blkmem, sdpisolver);
332 
333  return SCIP_OKAY;
334 }
335 
338  SCIP_SDPISOLVER* sdpisolver
339  )
340 {
341  assert( sdpisolver != NULL );
342 
343  sdpisolver->sdpcounter++;
344 
345  return SCIP_OKAY;
346 }
347 
350  SCIP_SDPISOLVER* sdpisolver
351  )
352 {
353  assert( sdpisolver != NULL );
354 
355  SCIPdebugMessage("Resetting counter of SDP-Interface from %d to 0.\n", sdpisolver->sdpcounter);
356  sdpisolver->sdpcounter = 0;
357 
358  return SCIP_OKAY;
359 }
360 
364 /*
365  * Solving Methods
366  */
367 
384  SCIP_SDPISOLVER* sdpisolver,
385  int nvars,
386  SCIP_Real* obj,
387  SCIP_Real* lb,
388  SCIP_Real* ub,
389  int nsdpblocks,
390  int* sdpblocksizes,
391  int* sdpnblockvars,
392  int sdpconstnnonz,
393  int* sdpconstnblocknonz,
395  int** sdpconstrow,
396  int** sdpconstcol,
397  SCIP_Real** sdpconstval,
398  int sdpnnonz,
399  int** sdpnblockvarnonz,
401  int** sdpvar,
403  int*** sdprow,
404  int*** sdpcol,
405  SCIP_Real*** sdpval,
406  int** indchanges,
408  int* nremovedinds,
409  int* blockindchanges,
410  int nremovedblocks,
411  int nlpcons,
412  int noldlpcons,
413  SCIP_Real* lplhs,
414  SCIP_Real* lprhs,
415  int* rownactivevars,
416  int lpnnonz,
417  int* lprow,
418  int* lpcol,
419  SCIP_Real* lpval,
420  SCIP_Real* start
421  )
422 {
423  return SCIPsdpiSolverLoadAndSolveWithPenalty(sdpisolver, 0.0, TRUE, TRUE, nvars, obj, lb, ub, nsdpblocks, sdpblocksizes, sdpnblockvars,
424  sdpconstnnonz, sdpconstnblocknonz, sdpconstrow, sdpconstcol, sdpconstval, sdpnnonz, sdpnblockvarnonz, sdpvar, sdprow, sdpcol, sdpval,
425  indchanges, nremovedinds, blockindchanges, nremovedblocks, nlpcons, noldlpcons, lplhs, lprhs, rownactivevars, lpnnonz, lprow, lpcol,
426  lpval, start, NULL);
427 }
428 
449  SCIP_SDPISOLVER* sdpisolver,
450  SCIP_Real penaltyparam,
451  SCIP_Bool withobj,
452  SCIP_Bool rbound,
453  int nvars,
454  SCIP_Real* obj,
455  SCIP_Real* lb,
456  SCIP_Real* ub,
457  int nsdpblocks,
458  int* sdpblocksizes,
459  int* sdpnblockvars,
460  int sdpconstnnonz,
461  int* sdpconstnblocknonz,
463  int** sdpconstrow,
464  int** sdpconstcol,
465  SCIP_Real** sdpconstval,
466  int sdpnnonz,
467  int** sdpnblockvarnonz,
469  int** sdpvar,
471  int*** sdprow,
472  int*** sdpcol,
473  SCIP_Real*** sdpval,
474  int** indchanges,
476  int* nremovedinds,
477  int* blockindchanges,
478  int nremovedblocks,
479  int nlpcons,
480  int noldlpcons,
481  SCIP_Real* lplhs,
482  SCIP_Real* lprhs,
483  int* rownactivevars,
484  int lpnnonz,
485  int* lprow,
486  int* lpcol,
487  SCIP_Real* lpval,
488  SCIP_Real* start,
489  SCIP_Bool* feasorig
491 )
492 {/*lint --e{413}*/
493  int* dsdpconstind = NULL; /* indices for constant SDP-constraint-matrices, needs to be stored for DSDP during solving and be freed only afterwards */
494  double* dsdpconstval = NULL; /* non-zero values for constant SDP-constraint-matrices, needs to be stored for DSDP during solving and be freed only afterwards */
495  int* dsdpind = NULL; /* indices for SDP-constraint-matrices, needs to be stored for DSDP during solving and be freed only afterwards */
496  double* dsdpval = NULL; /* non-zero values for SDP-constraint-matrices, needs to be stored for DSDP during solving and be freed only afterwards */
497  int* dsdplpbegcol = NULL; /* starting-indices for all columns in LP, needs to be stored for DSDP during solving and be freed only afterwards */
498  int* dsdplprow = NULL; /* row indices in LP, needs to be stored for DSDP during solving and be freed only afterwards */
499  double* dsdplpval = NULL; /* nonzeroes in LP, needs to be stored for DSDP during solving and be freed only afterwards */
500  int i;
501  int j;
502  int ind;
503  int block;
504  int startind;
505  int nfixedvars;
506  int dsdpnlpnonz = 0;
507  int nrnonz = 0;
508 
509 #ifdef SCIP_DEBUG
510  DSDPTerminationReason reason; /* this will later be used to check if DSDP converged */
511 #endif
512 
513  assert( sdpisolver != NULL );
514  assert( penaltyparam > -1 * sdpisolver->epsilon );
515  assert( penaltyparam < sdpisolver->epsilon || feasorig != NULL );
516  assert( nvars > 0 );
517  assert( obj != NULL );
518  assert( lb != NULL );
519  assert( ub != NULL );
520  assert( nsdpblocks >= 0 );
521  assert( nsdpblocks == 0 || sdpblocksizes != NULL );
522  assert( nsdpblocks == 0 || sdpnblockvars != NULL );
523  assert( sdpconstnnonz >= 0 );
524  assert( nsdpblocks == 0 || sdpconstnblocknonz != NULL );
525  assert( nsdpblocks == 0 || sdpconstrow != NULL );
526  assert( nsdpblocks == 0 || sdpconstcol != NULL );
527  assert( nsdpblocks == 0 || sdpconstval != NULL );
528  assert( sdpnnonz >= 0 );
529  assert( nsdpblocks == 0 || sdpnblockvarnonz != NULL );
530  assert( nsdpblocks == 0 || sdpvar != NULL );
531  assert( nsdpblocks == 0 || sdprow != NULL );
532  assert( nsdpblocks == 0 || sdpcol != NULL );
533  assert( nsdpblocks == 0 || sdpval != NULL );
534  assert( nsdpblocks == 0 || indchanges != NULL );
535  assert( nsdpblocks == 0 || nremovedinds != NULL );
536  assert( nsdpblocks == 0 || blockindchanges != NULL );
537  assert( 0 <= nremovedblocks && nremovedblocks <= nsdpblocks );
538  assert( nlpcons >= 0 );
539  assert( noldlpcons >= nlpcons );
540  assert( nlpcons == 0 || lplhs != NULL );
541  assert( nlpcons == 0 || lprhs != NULL );
542  assert( nlpcons == 0 || rownactivevars != NULL );
543  assert( lpnnonz >= 0 );
544  assert( nlpcons == 0 || lprow != NULL );
545  assert( nlpcons == 0 || lpcol != NULL );
546  assert( nlpcons == 0 || lpval != NULL );
547 
548  /* 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
549  * same SDP) */
550  if ( penaltyparam < sdpisolver->epsilon )
551  SCIPdebugMessage("Inserting Data into DSDP for SDP (%d) \n", ++sdpisolver->sdpcounter);
552  else
553  SCIPdebugMessage("Inserting Data again into DSDP for SDP (%d) \n", sdpisolver->sdpcounter);
554 
555  /* allocate memory for inputtodsdpmapper, dsdptoinputmapper and the fixed variable information, for the latter this will
556  * later be shrinked if the needed size is known */
557  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->inputtodsdpmapper), sdpisolver->nvars, nvars) );
558  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->dsdptoinputmapper), sdpisolver->nactivevars, nvars) );
559  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->fixedvarsval), sdpisolver->nvars - sdpisolver->nactivevars, nvars) ); /*lint !e776*/
560 
561  sdpisolver->nvars = nvars;
562  sdpisolver->nactivevars = 0;
563  nfixedvars = 0;
564 
565  /* find the fixed variables */
566  sdpisolver->fixedvarsobjcontr = 0.0;
567  for (i = 0; i < nvars; i++)
568  {
569  if ( isFixed(sdpisolver, lb[i], ub[i]) )
570  {
571  nfixedvars++;
572  sdpisolver->inputtodsdpmapper[i] = -nfixedvars;
573  sdpisolver->fixedvarsobjcontr += obj[i] * lb[i]; /* this is the value this variable contributes to the objective */
574  sdpisolver->fixedvarsval[nfixedvars - 1] = lb[i]; /* if lb=ub, than this is the value the variable will have in every solution */
575  SCIPdebugMessage("Fixing variable %d locally to %f for SDP %d in DSDP\n", i, lb[i], sdpisolver->sdpcounter);
576  }
577  else
578  {
579  sdpisolver->dsdptoinputmapper[sdpisolver->nactivevars] = i;
580  sdpisolver->nactivevars++;
581  sdpisolver->inputtodsdpmapper[i] = sdpisolver->nactivevars; /* dsdp starts counting at 1, so we do this after increasing nactivevars */
582  SCIPdebugMessage("Variable %d becomes variable %d for SDP %d in DSDP\n", i, sdpisolver->inputtodsdpmapper[i], sdpisolver->sdpcounter);
583  }
584  }
585  assert( sdpisolver->nactivevars + nfixedvars == sdpisolver->nvars );
586  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
587  {
588  SCIPdebugMessage("Variable %d is the slack variable for the explicit penalty formulation\n", sdpisolver->nactivevars + 1);
589  }
590 
591  /* if we want to solve without objective, we reset fixedvarsobjcontr */
592  if ( ! withobj )
593  sdpisolver->fixedvarsobjcontr = 0.0;
594 
595  /* shrink the fixedvars and dsdptoinputmapper arrays to the right size */
596  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->fixedvarsval), nvars, nfixedvars) );
597  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->dsdptoinputmapper), nvars, sdpisolver->nactivevars) );
598 
599  /* insert data */
600 
601  if ( sdpisolver->dsdp != NULL )
602  {
603  DSDP_CALL( DSDPDestroy(sdpisolver->dsdp) ); /* if there already exists a DSDP-instance, destroy the old one */
604  }
605 
606  /* in case we don't want to bound r, we can't use the penalty formulation in DSDP and have to give r explicitly */
607  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
608  {
609  DSDP_CALLM( DSDPCreate(sdpisolver->nactivevars + 1, &(sdpisolver->dsdp)) );
610  sdpisolver->penaltyworbound = TRUE;
611  }
612  else
613  {
614  DSDP_CALLM( DSDPCreate(sdpisolver->nactivevars, &(sdpisolver->dsdp)) );
615  sdpisolver->penaltyworbound = FALSE;
616  }
617  DSDP_CALLM( DSDPCreateSDPCone(sdpisolver->dsdp, nsdpblocks - nremovedblocks, &(sdpisolver->sdpcone)) );
618  DSDP_CALLM( DSDPCreateLPCone(sdpisolver->dsdp, &(sdpisolver->lpcone)) );
619  DSDP_CALLM( DSDPCreateBCone(sdpisolver->dsdp, &(sdpisolver->bcone)) );
620 
621 #ifdef SCIP_MORE_DEBUG
622  printf("setting objective values for SDP %d:\n", sdpisolver->sdpcounter);
623 #endif
624  for (i = 0; i < sdpisolver->nactivevars; i++)
625  {
626  if ( withobj )
627  {
628  /* insert objective value, DSDP counts from 1 to n instead of 0 to n-1, *(-1) because DSDP maximizes instead of minimizing */
629  DSDP_CALL( DSDPSetDualObjective(sdpisolver->dsdp, i+1, -1.0 * obj[sdpisolver->dsdptoinputmapper[i]]) );
630 #ifdef SCIP_MORE_DEBUG
631  printf("var %d (was var %d): %f, ", i+1, sdpisolver->dsdptoinputmapper[i], obj[sdpisolver->dsdptoinputmapper[i]]);
632 #endif
633  }
634  else
635  {
636  DSDP_CALL( DSDPSetDualObjective(sdpisolver->dsdp, i+1, 0.0) );
637  }
638  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, lb[sdpisolver->dsdptoinputmapper[i]]) )
639  {
640  /* insert lower bound, DSDP counts from 1 to n instead of 0 to n-1 */
641  DSDP_CALL( BConeSetLowerBound(sdpisolver->bcone, i+1, lb[sdpisolver->dsdptoinputmapper[i]]) );
642  }
643  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, ub[sdpisolver->dsdptoinputmapper[i]]) )
644  {
645  /* insert upper bound, DSDP counts from 1 to n instead of 0 to n-1 */
646  DSDP_CALL(BConeSetUpperBound(sdpisolver->bcone, i+1, ub[sdpisolver->dsdptoinputmapper[i]]));
647  }
648  }
649 
650  /* insert the objective value for r if solving without rbound, it is variable nactivevars + 1 and the objective is multiplied by -1 as we maximize */
651  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
652  {
653  DSDP_CALL( DSDPSetDualObjective(sdpisolver->dsdp, sdpisolver->nactivevars + 1, -1.0 * penaltyparam) );
654 #ifdef SCIP_MORE_DEBUG
655  printf("slack variable r: %f, ", penaltyparam);
656 #endif
657  }
658 
659 #ifdef SCIP_MORE_DEBUG
660  printf("\n");
661  SCIPdebugMessage("ATTENTION: BConeView shows the WRONG sign for the lower bound!\n");
662  BConeView(sdpisolver->bcone);
663 #endif
664 
665  /* set blocksizes */
666  for (block = 0; block < nsdpblocks; ++block)
667  {
668  /* only insert blocksizes for the blocks we didn't remove */
669  if ( blockindchanges[block] > -1 )
670  {
671  /* (blocks are counted from 0 to m-1) */
672  DSDP_CALL( SDPConeSetBlockSize(sdpisolver->sdpcone, block- blockindchanges[block], sdpblocksizes[block] - nremovedinds[block]) );
673  }
674  }
675 
676  /* start inserting the non-constant SDP-Constraint-Matrices */
677  if ( sdpnnonz > 0 )
678  {
679  int v;
680  int k;
681  int blockvar;
682 
683  nrnonz = 0;
684 
685  /* allocate memory */
686  /* This needs to be one long array, because DSDP uses it for solving so all nonzeros have to be in it and it may not be freed before the
687  * problem is solved. The distinct blocks/variables (for the i,j-parts) are then given by dsdpind + startind, which gives a pointer to the
688  * first array-element belonging to this block and then the number of elements in this block is given to DSDP for iterating over it. */
689 
690  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
691  {
692  /* we need to compute the total number of nonzeros for the slack variable r, which equals the total number of diagonal entries */
693  for (block = 0; block < nsdpblocks; block++)
694  {
695  nrnonz += sdpblocksizes[block] - nremovedinds[block];
696  }
697  assert( nrnonz >= 0 );
698 
699  /* indices given to DSDP, for this the elements in the lower triangular part of the matrix are labeled from 0 to n*(n+1)/2 -1 */
700  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdpind, sdpnnonz + nrnonz) ); /*lint !e776*/
701  /* values given to DSDP, these will be multiplied by -1 because in DSDP -1* (sum A_i^j y_i - A_0) should be positive semidefinite */
702  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdpval, sdpnnonz + nrnonz) ); /*lint !e776*/
703  }
704  else
705  {
706  /* indices given to DSDP, for this the elements in the lower triangular part of the matrix are labeled from 0 to n*(n+1)/2 -1 */
707  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdpind, sdpnnonz) );
708  /* values given to DSDP, these will be multiplied by -1 because in DSDP -1* (sum A_i^j y_i - A_0) should be positive semidefinite */
709  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdpval, sdpnnonz) );
710  }
711 
712  ind = 0; /* this will be used for iterating over the nonzeroes */
713 
714  for (block = 0; block < nsdpblocks; block++)
715  {
716  for (i = 0; i < sdpisolver->nactivevars; i++)
717  {
718  /* we iterate over all non-fixed variables, so add them to the dsdp arrays for this block/var combination */
719  v = sdpisolver->dsdptoinputmapper[i];
720 
721  /* find the position of variable v in this block */
722  blockvar = -1;
723  for (k = 0; k < sdpnblockvars[block]; k++)
724  {
725  if ( v == sdpvar[block][k] )
726  {
727  blockvar = k;
728  break;
729  }
730  }
731 
732  startind = ind;
733 
734  if ( blockvar > -1 ) /* the variable exists in this block */
735  {
736  for (k = 0; k < sdpnblockvarnonz[block][blockvar]; k++)
737  {
738  /* rows and cols with active nonzeros should not be removed */
739  assert( indchanges[block][sdprow[block][blockvar][k]] > -1 && indchanges[block][sdpcol[block][blockvar][k]] > -1 );
740 
741  /* substract the number of removed indices before the row and col to get the indices after fixings */
742  dsdpind[ind] = compLowerTriangPos(sdprow[block][blockvar][k] - indchanges[block][sdprow[block][blockvar][k]],
743  sdpcol[block][blockvar][k] - indchanges[block][sdpcol[block][blockvar][k]]);
744  dsdpval[ind] = -1.0 * sdpval[block][blockvar][k]; /* *(-1) because in DSDP -1* (sum A_i^j y_i - A_0) should be positive semidefinite */
745  ind++;
746  }
747 
748  /* sort the arrays for this Matrix (by non decreasing indices) as this might help the solving time of DSDP */
749  SCIPsortIntReal(dsdpind + startind, dsdpval + startind, sdpnblockvarnonz[block][blockvar]);
750 
751  assert( blockindchanges[block] > -1 ); /* we shouldn't insert into blocks we removed */
752 
753  /* i + 1 because DSDP starts counting the variables at 1, adding startind shifts the arrays to the first
754  * nonzero belonging to this block and this variable */
755  DSDP_CALL( SDPConeSetASparseVecMat(sdpisolver->sdpcone, block - blockindchanges[block], i + 1, sdpblocksizes[block] - nremovedinds[block],
756  1.0, 0, dsdpind + startind,dsdpval + startind, sdpnblockvarnonz[block][blockvar]));
757  }
758  }
759  }
760  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
761  {
762  startind = ind;
763  /* add r * Identity for each block */
764  for (block = 0; block < nsdpblocks; block++)
765  {
766  if ( blockindchanges[block] > -1 )
767  {
768  for (i = 0; i < sdpblocksizes[block] - nremovedinds[block]; i++)
769  {
770  dsdpind[ind] = compLowerTriangPos(i, i);
771  dsdpval[ind] = -1.0; /* *(-1) because in DSDP -1* (sum A_i^j y_i - A_0 + r*I) should be positive semidefinite */
772  ind++;
773  }
774  DSDP_CALL( SDPConeSetASparseVecMat(sdpisolver->sdpcone, block - blockindchanges[block], sdpisolver->nactivevars + 1,
775  sdpblocksizes[block] - nremovedinds[block], 1.0, 0, dsdpind + ind - (sdpblocksizes[block] - nremovedinds[block]) ,
776  dsdpval + ind - (sdpblocksizes[block] - nremovedinds[block]), sdpblocksizes[block] - nremovedinds[block]) ); /*lint !e679*/
777  }
778  }
779  assert( ind - startind == nrnonz );
780  }
781  }
782 
783  /* start inserting the constant matrix */
784  if ( sdpconstnnonz > 0 )
785  {
786  assert( nsdpblocks > 0 );
787  assert( sdpconstnblocknonz!= NULL );
788  assert( sdpconstcol != NULL );
789  assert( sdpconstrow != NULL );
790  assert( sdpconstval != NULL );
791 
792  /* allocate memory */
793 
794  /* DSDP uses these for solving, so they may not be freed before the problem is solved. */
795 
796  /* indices given to DSDP, for this the elements in the lower triangular part of the matrix are labeled from 0 to n*(n+1)/2 -1 */
797  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdpconstind, sdpconstnnonz) );
798  /* values given to DSDP, for this the original values are mutliplied by -1 because in DSDP -1* (sum A_i^j y_i - A_0) should be positive semidefinite */
799  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdpconstval, sdpconstnnonz) );
800 
801  ind = 0;
802 
803  for (block = 0; block < nsdpblocks; block++)
804  {
805  startind = ind; /* starting index of this block in the dsdpconst arrays */
806 
807  if ( sdpconstnblocknonz[block] > 0 )
808  {
809  /* insert the constant-nonzeros */
810  for (i = 0; i < sdpconstnblocknonz[block]; i++)
811  {
812  /* rows and cols with nonzeros should not be removed */
813  assert( indchanges[block][sdpconstrow[block][i]] > -1 && indchanges[block][sdpconstcol[block][i]] > -1 );
814 
815  /* substract the number of deleted indices before this to get the index after variable fixings */
816  dsdpconstind[ind] = compLowerTriangPos(sdpconstrow[block][i] - indchanges[block][sdpconstrow[block][i]],
817  sdpconstcol[block][i] - indchanges[block][sdpconstcol[block][i]]);
818  dsdpconstval[ind] = -1 * sdpconstval[block][i]; /* *(-1) because in DSDP -1* (sum A_i^j y_i - A_0^j) should be positive semidefinite */
819  ind++;
820  }
821 
822  /* sort the arrays for this Matrix (by non decreasing indices) as this might help the solving time of DSDP */
823  SCIPsortIntReal(dsdpconstind + startind, dsdpconstval + startind, sdpconstnblocknonz[block]);
824 
825  assert( blockindchanges[block] > -1 ); /* we shouldn't insert into a block we removed */
826 
827  /* constant matrix is given as variable 0, the arrays are shifted to the first element of this block by adding
828  * startind, ind - startind gives the number of elements for this block */
829  DSDP_CALL( SDPConeSetASparseVecMat(sdpisolver->sdpcone, block - blockindchanges[block], 0, sdpblocksizes[block] - nremovedinds[block],
830  1.0, 0, dsdpconstind + startind, dsdpconstval + startind, ind - startind));
831  }
832  }
833  }
834 
835 #ifdef SCIP_MORE_DEBUG
836  SDPConeView2(sdpisolver->sdpcone);
837 #endif
838 
839  /* start inserting the LP constraints */
840  if ( nlpcons > 0 || lpnnonz > 0 || ! SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit) )
841  {
842  int nextcol;
843  int* rowmapper; /* maps the lhs- and rhs-inequalities of the old LP-cons to their constraint numbers in DSDP */
844  int pos;
845  int newpos;
846  int nlpineqs;
847 
848  assert( noldlpcons > 0 );
849  assert( lprhs != NULL );
850  assert( lpcol != NULL );
851  assert( lprow != NULL );
852  assert( lpval != NULL );
853 
854  /* 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 */
855  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &rowmapper, 2*noldlpcons) ); /*lint !e647*/
856 
857  /* compute the rowmapper and the number of inequalities we have to add to DSDP (as we have to split the ranged rows) */
858  pos = 0;
859  newpos = 0; /* the position in the lhs and rhs arrays */
860  for (i = 0; i < noldlpcons; i++)
861  {
862  if ( rownactivevars[i] >= 2 )
863  {
864  if ( lplhs[newpos] > - SCIPsdpiSolverInfinity(sdpisolver) )
865  {
866  rowmapper[2*i] = pos; /*lint !e679*/
867  pos++;
868  }
869  else
870  rowmapper[2*i] = -1; /*lint !e679*/
871  if ( lprhs[newpos] < SCIPsdpiSolverInfinity(sdpisolver) )
872  {
873  rowmapper[2*i + 1] = pos; /*lint !e679*/
874  pos++;
875  }
876  else
877  rowmapper[2*i + 1] = -1; /*lint !e679*/
878 
879  newpos++;
880  }
881  else
882  {
883  rowmapper[2*i] = -1; /*lint !e679*/
884  rowmapper[2*i + 1] = -1; /*lint !e679*/
885  }
886  }
887  nlpineqs = pos;
888  assert( nlpineqs <= 2*nlpcons ); /* *2 comes from left- and right-hand-sides */
889 
890  /* memory allocation */
891 
892  /* these arrays are needed in DSDP during solving, so they may only be freed afterwards */
893  /* dsdplpbegcol[i] gives the number of nonzeros in column 0 (right hand side) till i-1 (i going from 1 till m, with extra entries 0 (always 0)
894  * and m+1 (always lpcons + lpnnonz)) */
895  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
896  {
897  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdplpbegcol, sdpisolver->nactivevars + 3) ); /* extra entry for r */ /*lint !e776*/
898  }
899  else
900  {
901  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdplpbegcol, sdpisolver->nactivevars + 2) ); /*lint !e776*/
902  }
903  /* dsdplprow saves the row indices of the LP for DSDP */
904  /* worst-case length is 2*lpnnonz + nlpineqs, because left- and right-hand-sides are also included in the vectorand we might have to duplicate the
905  * non-zeros when splitting the ranged rows, this will be shortened after the exact length after fixings is known, in case we have an objective limit,
906  * this is increased by one entry for the right-hand-side and at most nvars entries for the nonzeros, in case of rbound = FALSE, where we have to add
907  * the entries for r ourselves, we have to add another nlpineqs for one entry for r for each active lp-constraint */
908  if ( SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit) )
909  {
910  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
911  {
912  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdplprow, 2 * nlpineqs + 2*lpnnonz) ); /*lint !e647*/
913  }
914  else
915  {
916  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdplprow, nlpineqs + 2*lpnnonz) ); /*lint !e647*/
917  }
918  }
919  else
920  {
921  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
922  {
923  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdplprow, (nlpineqs + 1) + 2*lpnnonz + nvars + nlpineqs) ); /*lint !e647*/
924  }
925  else
926  {
927  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdplprow, (nlpineqs + 1) + 2*lpnnonz + nvars) ); /*lint !e647*/
928  }
929  }
930 
931  /* values given to DSDP */
932  /* dsdplprow saves the row indices of the LP for DSDP */
933  /* worst-case length is 2*lpnnonz + nlpineqs, because left- and right-hand-sides are also included in the vectorand we might have to duplicate the
934  * non-zeros when splitting the ranged rows, this will be shortened after the exact length after fixings is known, in case we have an objective limit,
935  * this is increased by one entry for the right-hand-side and at most nvars entries for the nonzeros, in case of rbound = FALSE, where we have to add
936  * the entries for r ourselves, we have to add another nlpineqs for one entry for r for each active lp-constraint */
937  if ( SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit) )
938  {
939  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
940  {
941  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdplpval, 2 * nlpineqs + 2*lpnnonz) ); /*lint !e647*/
942  }
943  else
944  {
945  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdplpval, nlpineqs + 2*lpnnonz) ); /*lint !e647*/
946  }
947  }
948  else
949  {
950  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
951  {
952  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdplpval, (nlpineqs + 1) + 2*lpnnonz + nvars + nlpineqs) ); /*lint !e647*/
953  }
954  else
955  {
956  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdplpval, (nlpineqs + 1) + 2*lpnnonz + nvars) ); /*lint !e647*/
957  }
958  }
959 
960  /* add all left- and right-hand-sides that are greater than zero (if their corresponding inequalities exist), the pos counter is increased for every
961  * active row, to get the correct row numbers, the nonz-counter only if the lhs/rhs is unequal to zero and added to DSDP */
962  dsdpnlpnonz = 0;
963  pos = 0;
964  for (i = 0; i < nlpcons; i++)
965  {
966  if ( lplhs[i] > - SCIPsdpiSolverInfinity(sdpisolver) )
967  {
968  if ( REALABS(lplhs[i]) > sdpisolver->epsilon )
969  {
970  dsdplprow[dsdpnlpnonz] = pos;
971  dsdplpval[dsdpnlpnonz] = -lplhs[i]; /* we multiply by -1 because DSDP wants <= instead of >= */
972  dsdpnlpnonz++;
973  }
974  pos++;
975  }
976  if ( lprhs[i] < SCIPsdpiSolverInfinity(sdpisolver) )
977  {
978  if ( REALABS(lprhs[i]) > sdpisolver->epsilon )
979  {
980  dsdplprow[dsdpnlpnonz] = pos;
981  dsdplpval[dsdpnlpnonz] = lprhs[i];
982  dsdpnlpnonz++;
983  }
984  pos++;
985  }
986  }
987  assert( pos == nlpineqs );
988 
989  /* add the right-hand-side for the objective bound */
990  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit) )
991  {
992  if ( REALABS(sdpisolver->objlimit) > sdpisolver->epsilon )
993  {
994  dsdplprow[dsdpnlpnonz] = nlpcons; /* this is the last lp-constraint, as DSDP counts from 0 to nlpcons-1, this is number nlpcons */
995  dsdplpval[dsdpnlpnonz] = sdpisolver->objlimit; /* as we want <= upper bound, this is the correct type of inequality for DSDP */
996  dsdpnlpnonz++;
997  }
998  }
999 
1000  /* now add the nonzeros */
1001 
1002  /* for this we have to sort the nonzeros by col first and then by row, as this is the sorting DSDP wants */
1003  sortColRow(lprow, lpcol, lpval, lpnnonz);
1004 
1005  /* iterate over all nonzeros to add the active ones to the dsdp arrays and compute dsdplpbegcol */
1006  nextcol = 0;
1007  dsdplpbegcol[0] = 0;
1008  for (i = 0; i < lpnnonz; i++)
1009  {
1010  /* if a new variable starts, set the corresponding dsdplpbegcol-entry */
1011  if ( lpcol[i] >= nextcol )
1012  {
1013  /* set the dsdplpbegcol entries, as there might be active variables which appear only in the sdp but not the lp-part, we also have to set
1014  * the starting values for all variables in between to the same value (as we also set the entry for the found variable, this for-queue
1015  * will always have at least one index in the index set) */
1016  for (j = nextcol; j <= lpcol[i]; j++)
1017  {
1018  if ( sdpisolver->inputtodsdpmapper[j] >= 0 )
1019  {
1020  assert( ! (isFixed(sdpisolver, lb[j], ub[j])) );
1021  dsdplpbegcol[sdpisolver->inputtodsdpmapper[j]] = dsdpnlpnonz;
1022 
1023  /* add the entry to the objlimit-lp-constraint for the last variables */
1024  if ( (! SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit)) && (REALABS( obj[j] ) > sdpisolver->epsilon))
1025  {
1026  dsdplprow[dsdpnlpnonz] = nlpcons;
1027  dsdplpval[dsdpnlpnonz] = obj[j];
1028  dsdpnlpnonz++;
1029  }
1030  }
1031  }
1032  nextcol = j;
1033  }
1034  /* add the nonzero, if it isn't fixed and the row isn't to be deleted (because it is only a variable bound) */
1035  if ( ! isFixed(sdpisolver, lb[lpcol[i]], ub[lpcol[i]]) )
1036  {
1037  /* add it to the inequality for the lhs of the ranged row, if it exists */
1038  if ( rowmapper[2*lprow[i]] > -1 ) /*lint !e679*/
1039  {
1040  /* the index is adjusted for deleted lp rows, also rows are numbered 0,...,nlpcons-1 in DSDP, as they are
1041  * here, nlpcons is added to the index as the first nlpcons entries correspond to the right hand sides */
1042  dsdplprow[dsdpnlpnonz] = rowmapper[2*lprow[i]]; /*lint !e679*/
1043  dsdplpval[dsdpnlpnonz] = -lpval[i]; /* - because dsdp wants <= instead of >= constraints */
1044  dsdpnlpnonz++;
1045  }
1046  /* add it to the inequality for the rhs of the ranged row, if it exists */
1047  if ( rowmapper[2*lprow[i] + 1] > -1 ) /*lint !e679*/
1048  {
1049  /* the index is adjusted for deleted lp rows, also rows are numbered 0,...,nlpcons-1 in DSDP, as they are
1050  * here, nlpcons is added to the index as the first nlpcons entries correspond to the right hand sides */
1051  dsdplprow[dsdpnlpnonz] = rowmapper[2*lprow[i] + 1]; /*lint !e679*/
1052  dsdplpval[dsdpnlpnonz] = lpval[i];
1053  dsdpnlpnonz++;
1054  }
1055  }
1056 #ifndef SCIP_NDEBUG
1057  /* if this is an active nonzero for the row, it should have at least one active var */
1058  else
1059  assert( isFixed(sdpisolver, lb[lpcol[i]], ub[lpcol[i]]) || rownactivevars[lprow[i]] == 1 );
1060 #endif
1061  }
1062 
1063  /* set the begcol array for all remaining variables (that aren't part of the LP-part), and also set the objlimit-constraint-entries */
1064  for (j = nextcol; j < nvars; j++)
1065  {
1066  if ( sdpisolver->inputtodsdpmapper[j] >= 0 )
1067  {
1068  assert( ! (isFixed(sdpisolver, lb[j], ub[j])) );
1069  dsdplpbegcol[sdpisolver->inputtodsdpmapper[j]] = dsdpnlpnonz;
1070  /* add the entry to the objlimit-lp-constraint for the last variables */
1071  if ( (! SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit)) && (REALABS( obj[j] ) > sdpisolver->epsilon))
1072  {
1073  dsdplprow[dsdpnlpnonz] = nlpcons;
1074  dsdplpval[dsdpnlpnonz] = obj[j];
1075  dsdpnlpnonz++;
1076  }
1077  }
1078  }
1079 
1080  dsdplpbegcol[sdpisolver->nactivevars + 1] = dsdpnlpnonz; /*lint !e679*/
1081 
1082  /* add r * Identity if using a penalty formulation without a bound on r */
1083  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
1084  {
1085  for (i = 0; i < nlpineqs; i++)
1086  {
1087  dsdplprow[dsdpnlpnonz] = i;
1088  dsdplpval[dsdpnlpnonz] = -1.0; /* for >=-inequalities we would add a +1, but then we have to multiply these with -1 for DSDP */
1089  dsdpnlpnonz++;
1090  }
1091  dsdplpbegcol[sdpisolver->nactivevars + 2] = dsdpnlpnonz; /*lint !e679*/
1092  }
1093 
1094  /* free the memory for the rowshifts-array */
1095  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &rowmapper, 2*noldlpcons); /*lint !e647*/
1096 
1097  /* shrink the dsdplp-arrays */
1098  if ( SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit) )
1099  {
1100  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
1101  {
1102  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &dsdplprow, 2*nlpineqs + 2*lpnnonz, dsdpnlpnonz) ); /*lint !e647*/
1103  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &dsdplpval, 2*nlpineqs + 2*lpnnonz, dsdpnlpnonz) ); /*lint !e647*/
1104  }
1105  else
1106  {
1107  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &dsdplprow, nlpineqs + 2*lpnnonz, dsdpnlpnonz) ); /*lint !e647*/
1108  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &dsdplpval, nlpineqs + 2*lpnnonz, dsdpnlpnonz) ); /*lint !e647*/
1109  }
1110  }
1111  else
1112  {
1113  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
1114  {
1115  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &dsdplprow, (nlpineqs + 1) + 2*lpnnonz + nvars + nlpineqs, dsdpnlpnonz) ); /*lint !e647*/
1116  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &dsdplpval, (nlpineqs + 1) + 2*lpnnonz + nvars + nlpineqs, dsdpnlpnonz) ); /*lint !e647*/
1117  }
1118  else
1119  {
1120  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &dsdplprow, (nlpineqs + 1) + 2*lpnnonz + nvars, dsdpnlpnonz) ); /*lint !e647*/
1121  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &dsdplpval, (nlpineqs + 1) + 2*lpnnonz + nvars, dsdpnlpnonz) ); /*lint !e647*/
1122  }
1123  }
1124  /* add the arrays to dsdp (in this case we need no additional if for the penalty version without bounds, as we already added the extra var,
1125  * so DSDP knows, that there is an additional entry in dsdplpbegcol which then gives the higher number of nonzeros) */
1126  if ( SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit) )
1127  {
1128  DSDP_CALL( LPConeSetData(sdpisolver->lpcone, nlpineqs, dsdplpbegcol, dsdplprow, dsdplpval) );
1129  }
1130  else
1131  {
1132  DSDP_CALL( LPConeSetData(sdpisolver->lpcone, nlpineqs + 1, dsdplpbegcol, dsdplprow, dsdplpval) );
1133  }
1134 #ifdef SCIP_MORE_DEBUG
1135  LPConeView(sdpisolver->lpcone);
1136 #endif
1137  }
1138 
1139  SCIPdebugMessage("Calling DSDP-Solve for SDP (%d) \n", sdpisolver->sdpcounter);
1140 
1141  DSDP_CALL( DSDPSetGapTolerance(sdpisolver->dsdp, sdpisolver->epsilon) ); /* set DSDP's tolerance for duality gap */
1142  DSDP_CALL( DSDPSetRTolerance(sdpisolver->dsdp, sdpisolver->feastol) ); /* set DSDP's tolerance for the SDP-constraints */
1143  if ( sdpisolver-> sdpinfo )
1144  {
1145  DSDP_CALL( DSDPSetStandardMonitor(sdpisolver->dsdp, 1) ); /* output DSDP information after every 1 iteration */
1146  }
1147 
1148  /* set the penalty parameter (only if rbound = TRUE, otherwise we had to add everything ourselves) */
1149  if ( penaltyparam >= sdpisolver->epsilon && rbound ) /* in sdpisolverSolve this is called with an exact 0 */
1150  {
1151  DSDP_CALL( DSDPSetPenaltyParameter(sdpisolver->dsdp, penaltyparam) );
1152  DSDP_CALL( DSDPUsePenalty(sdpisolver->dsdp, 1) );
1153  }
1154 
1155  /* set the starting solution */
1156  if ( start != NULL )
1157  {
1158  for (i = 1; i <= sdpisolver->nactivevars; i++) /* we iterate over the variables in DSDP */
1159  {
1160  DSDP_CALL( DSDPSetY0(sdpisolver->dsdp, i, start[sdpisolver->dsdptoinputmapper[i]]) );
1161  }
1162  }
1163 
1164  /* start the solving process */
1165  DSDP_CALLM( DSDPSetup(sdpisolver->dsdp) );
1166  DSDP_CALL( DSDPSolve(sdpisolver->dsdp) );
1167 
1168  DSDP_CALL( DSDPComputeX(sdpisolver->dsdp) ); /* computes X and determines feasibility and unboundedness of the solution */
1169  sdpisolver->solved = TRUE;
1170 
1171  /*these arrays were used to give information to DSDP and were needed during solving and for computing X, so they may only be freed now*/
1172  if ( sdpconstnnonz > 0 )
1173  {
1174  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &dsdpconstval, sdpconstnnonz);
1175  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &dsdpconstind, sdpconstnnonz);
1176  }
1177 
1178  if ( sdpnnonz > 0 )
1179  {
1180  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
1181  {
1182  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &dsdpval, sdpnnonz + nrnonz); /*lint !e776*/
1183  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &dsdpind, sdpnnonz + nrnonz); /*lint !e776*/
1184  }
1185  else
1186  {
1187  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &dsdpval, sdpnnonz);
1188  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &dsdpind, sdpnnonz);
1189  }
1190  }
1191 
1192  if ( nlpcons > 0 || lpnnonz > 0 )
1193  {
1194  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &dsdplpval, dsdpnlpnonz); /*lint !e644*/
1195  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &dsdplprow, dsdpnlpnonz);
1196  if ( penaltyparam > sdpisolver->epsilon && (! rbound) )
1197  {
1198  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &dsdplpbegcol, sdpisolver->nactivevars + 3); /*lint !e776*/
1199  }
1200  else
1201  {
1202  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &dsdplpbegcol, sdpisolver->nactivevars + 2); /*lint !e776*/
1203  }
1204  }
1205 
1206 #ifdef SCIP_DEBUG
1207  DSDP_CALL(DSDPStopReason(sdpisolver->dsdp, &reason));
1208 
1209  switch ( reason )
1210  {
1211  case DSDP_CONVERGED:
1212  SCIPdebugMessage("DSDP converged!\n");
1213  break;
1214 
1215  case DSDP_INFEASIBLE_START:
1216  SCIPdebugMessage("DSDP started with an infeasible point!\n");
1217  break;
1218 
1219  case DSDP_SMALL_STEPS:
1220  SCIPdebugMessage("Short step lengths created by numerical difficulties prevented progress in DSDP!\n");
1221  break;
1222 
1223  case DSDP_INDEFINITE_SCHUR_MATRIX:
1224  SCIPdebugMessage("Schur Matrix in DSDP was indefinite but should have been positive semidefinite!\n");
1225  break;
1226 
1227  case DSDP_MAX_IT:
1228  SCIPdebugMessage("DSDP reached maximum number of iterations!\n");
1229  break;
1230 
1231  case DSDP_NUMERICAL_ERROR:
1232  SCIPdebugMessage("A numerical error occured in DSDP!\n");
1233  break;
1234 
1235  case DSDP_UPPERBOUND:
1236  SCIPdebugMessage("Dual objective value in DSDP reached upper bound.\n");
1237  break;
1238 
1239  case DSDP_USER_TERMINATION:
1240  SCIPdebugMessage("DSDP didn't stop solving, did you?\n");
1241  break;
1242 
1243  case CONTINUE_ITERATING:
1244  SCIPdebugMessage("DSDP wants to continue iterating but somehow was stopped!\n");
1245  break;
1246 
1247  default:
1248  SCIPdebugMessage("Unknown stopping reason in DSDP!\n");
1249  break;
1250  }
1251 #endif
1252 
1253  if ( penaltyparam >= sdpisolver->epsilon )
1254  {
1255  if ( rbound )
1256  {
1257  /* in this case we used the penalty-formulation of DSDP, so we can check their value of r */
1258  double rval;
1259 
1260  DSDP_CALL( DSDPGetR(sdpisolver->dsdp, &rval) );
1261 
1262  *feasorig = (rval < sdpisolver->feastol );
1263  if ( ! *feasorig )
1264  {
1265  SCIPdebugMessage("Solution not feasible in original problem, r = %f\n", rval);
1266  }
1267  }
1268  else
1269  {
1270  double* dsdpsol;
1271 
1272  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdpsol, sdpisolver->nactivevars + 1) ); /*lint !e776*/
1273  /* last entry of DSDPGetY needs to be the number of variables, will return an error otherwise */
1274  DSDP_CALL( DSDPGetY(sdpisolver->dsdp, dsdpsol, sdpisolver->nactivevars + 1) );
1275 
1276  *feasorig = (dsdpsol[sdpisolver->nactivevars] < sdpisolver->feastol); /* r is the last variable in DSDP, so the last entry gives us the value */
1277  if ( ! *feasorig )
1278  {
1279  SCIPdebugMessage("Solution not feasible in original problem, r = %f\n", dsdpsol[sdpisolver->nactivevars]);
1280  }
1281 
1282  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &dsdpsol, sdpisolver->nactivevars + 1); /*lint !e776*/
1283  }
1284  }
1285 
1286  return SCIP_OKAY;
1287 }
1293 /*
1294  * Solution Information Methods
1295  */
1296 
1302  SCIP_SDPISOLVER* sdpisolver
1303  )
1304 {
1305  assert( sdpisolver != NULL );
1306  return sdpisolver->solved;
1307 }
1308 
1316  SCIP_SDPISOLVER* sdpisolver
1317  )
1318 {
1319  DSDPSolutionType pdfeasible;
1320 
1321  assert( sdpisolver != NULL );
1322  CHECK_IF_SOLVED_BOOL( sdpisolver );
1323 
1324  DSDP_CALL_BOOL( DSDPGetSolutionType(sdpisolver->dsdp, &pdfeasible) );
1325 
1326  if ( pdfeasible == DSDP_PDUNKNOWN )
1327  return FALSE;
1328 
1329  return TRUE;
1330 }
1331 
1334  SCIP_SDPISOLVER* sdpisolver,
1335  SCIP_Bool* primalfeasible,
1336  SCIP_Bool* dualfeasible
1337  )
1338 {
1339  DSDPSolutionType pdfeasible;
1340 
1341  assert( sdpisolver != NULL );
1342  assert( primalfeasible != NULL );
1343  assert( dualfeasible != NULL );
1344  CHECK_IF_SOLVED( sdpisolver );
1345 
1346  DSDP_CALL( DSDPGetSolutionType(sdpisolver->dsdp, &pdfeasible) );
1347 
1348  switch ( pdfeasible )
1349  {
1350  case DSDP_PDFEASIBLE:
1351  *primalfeasible = TRUE;
1352  *dualfeasible = TRUE;
1353  break;
1354 
1355  case DSDP_UNBOUNDED:
1356  *primalfeasible = FALSE;
1357  *dualfeasible = TRUE;
1358  break;
1359 
1360  case DSDP_INFEASIBLE:
1361  *primalfeasible = TRUE;
1362  *dualfeasible = FALSE;
1363  break;
1364 
1365  default: /* should only include DSDP_PDUNKNOWN */
1366  SCIPerrorMessage("DSDP doesn't know if primal and dual solutions are feasible\n");
1367  return SCIP_LPERROR;
1368  }
1369 
1370  return SCIP_OKAY;
1371 }
1372 
1376  SCIP_SDPISOLVER* sdpisolver
1377  )
1378 {
1379  DSDPSolutionType pdfeasible;
1380 
1381  assert( sdpisolver != NULL );
1382  CHECK_IF_SOLVED_BOOL( sdpisolver );
1383 
1384  DSDP_CALL_BOOL( DSDPGetSolutionType(sdpisolver->dsdp, &pdfeasible) );
1385  if ( pdfeasible == DSDP_PDUNKNOWN )
1386  {
1387 /* SCIPerrorMessage("DSDP doesn't know if primal and dual solutions are feasible");
1388  SCIPABORT();
1389  return SCIP_LPERROR;*/
1390  SCIPdebugMessage("DSDP doesn't know if primal and dual solutions are feasible.");
1391  return FALSE;
1392  }
1393  else if ( pdfeasible == DSDP_INFEASIBLE )
1394  return TRUE;
1395 
1396  return FALSE;
1397 }
1398 
1402  SCIP_SDPISOLVER* sdpisolver
1403  )
1404 {
1405  DSDPSolutionType pdfeasible;
1406 
1407  assert( sdpisolver != NULL );
1408  CHECK_IF_SOLVED_BOOL( sdpisolver );
1409 
1410  DSDP_CALL_BOOL( DSDPGetSolutionType(sdpisolver->dsdp, &pdfeasible) );
1411  if ( pdfeasible == DSDP_PDUNKNOWN )
1412  {
1413 /* SCIPerrorMessage("DSDP doesn't know if primal and dual solutions are feasible");
1414  SCIPABORT();
1415  return SCIP_LPERROR;*/
1416  SCIPdebugMessage("DSDP doesn't know if primal and dual solutions are feasible");
1417  return FALSE;
1418  }
1419  else if ( pdfeasible == DSDP_UNBOUNDED )
1420  return TRUE;
1421 
1422  return FALSE;
1423 }
1424 
1428  SCIP_SDPISOLVER* sdpisolver
1429  )
1430 {
1431  DSDPSolutionType pdfeasible;
1432 
1433  assert( sdpisolver != NULL );
1434  CHECK_IF_SOLVED_BOOL( sdpisolver );
1435 
1436  DSDP_CALL_BOOL( DSDPGetSolutionType(sdpisolver->dsdp, &pdfeasible) );
1437  if ( pdfeasible == DSDP_PDUNKNOWN )
1438  {
1439  SCIPdebugMessage("DSDP doesn't know if primal and dual solutions are feasible");
1440  return FALSE;
1441  }
1442  else if ( pdfeasible == DSDP_UNBOUNDED )
1443  return FALSE;
1444 
1445  return TRUE;
1446 }
1447 
1451  SCIP_SDPISOLVER* sdpisolver
1452  )
1453 {
1454  DSDPSolutionType pdfeasible;
1455 
1456  assert( sdpisolver != NULL );
1457  CHECK_IF_SOLVED_BOOL( sdpisolver );
1458 
1459  DSDP_CALL_BOOL( DSDPGetSolutionType(sdpisolver->dsdp, &pdfeasible) );
1460  if ( pdfeasible == DSDP_PDUNKNOWN )
1461  {
1462  SCIPdebugMessage("DSDP doesn't know if primal and dual solutions are feasible");
1463  return FALSE;
1464  }
1465  else if ( pdfeasible == DSDP_UNBOUNDED )
1466  return TRUE;
1467 
1468  return FALSE;
1469 }
1470 
1474  SCIP_SDPISOLVER* sdpisolver
1475  )
1476 {
1477  DSDPSolutionType pdfeasible;
1478 
1479  assert( sdpisolver != NULL );
1480  CHECK_IF_SOLVED_BOOL( sdpisolver );
1481 
1482  DSDP_CALL_BOOL(DSDPGetSolutionType(sdpisolver->dsdp, &pdfeasible));
1483 
1484  if ( pdfeasible == DSDP_PDUNKNOWN )
1485  {
1486  SCIPdebugMessage("DSDP doesn't know if primal and dual solutions are feasible");
1487  return FALSE;
1488  }
1489  else if ( pdfeasible == DSDP_INFEASIBLE )
1490  return TRUE;
1491 
1492  return FALSE;
1493 }
1494 
1498  SCIP_SDPISOLVER* sdpisolver
1499  )
1500 {
1501  DSDPSolutionType pdfeasible;
1502 
1503  assert( sdpisolver != NULL );
1504  CHECK_IF_SOLVED_BOOL( sdpisolver );
1505 
1506  DSDP_CALL_BOOL( DSDPGetSolutionType(sdpisolver->dsdp, &pdfeasible) );
1507 
1508  if ( pdfeasible == DSDP_PDUNKNOWN )
1509  {
1510  SCIPdebugMessage("DSDP doesn't know if primal and dual solutions are feasible");
1511  return FALSE;
1512  }
1513  else if ( pdfeasible == DSDP_INFEASIBLE )
1514  return FALSE;
1515 
1516  return TRUE;
1517 }
1518 
1521  SCIP_SDPISOLVER* sdpisolver
1522  )
1523 {
1524  DSDPTerminationReason reason;
1525 
1526  assert( sdpisolver != NULL );
1527  CHECK_IF_SOLVED_BOOL( sdpisolver );
1528 
1529  DSDP_CALL_BOOL( DSDPStopReason(sdpisolver->dsdp, &reason) );
1530 
1531  if ( reason == DSDP_CONVERGED )
1532  return TRUE;
1533 
1534  return FALSE;
1535 }
1536 
1539  SCIP_SDPISOLVER* sdpisolver
1540  )
1541 {/*lint --e{715}*/
1542  SCIPdebugMessage("Method not implemented for DSDP, as objective limit is given as an ordinary LP-constraint, so in case the objective limit was "
1543  "exceeded, the problem will be reported as infeasible ! \n");
1544 
1545  return FALSE;
1546 }
1547 
1550  SCIP_SDPISOLVER* sdpisolver
1551  )
1552 {
1553  DSDPTerminationReason reason;
1554 
1555  assert( sdpisolver != NULL );
1556  CHECK_IF_SOLVED_BOOL( sdpisolver );
1557 
1558  DSDP_CALL_BOOL(DSDPStopReason(sdpisolver->dsdp, &reason));
1559 
1560  if ( reason == DSDP_MAX_IT )
1561  return TRUE;
1562 
1563  return FALSE;
1564 }
1565 
1568  SCIP_SDPISOLVER* sdpisolver
1569  )
1570 {/*lint --e{715}*/
1571  SCIPdebugMessage("Not implemented in DSDP!\n");
1572  return SCIP_LPERROR;
1573 }
1574 
1586  SCIP_SDPISOLVER* sdpisolver
1587  )
1588 {
1589  DSDPTerminationReason reason;
1590  int dsdpreturn;
1591 
1592  assert( sdpisolver != NULL );
1593 
1594  if ( sdpisolver->dsdp == NULL || (! sdpisolver->solved) )
1595  return -1;
1596 
1597  dsdpreturn = DSDPStopReason(sdpisolver->dsdp, &reason);
1598 
1599  if (dsdpreturn != 0)
1600  {
1601  SCIPerrorMessage("DSDP-Error <%d> in function call.\n", dsdpreturn);
1602  return 7;
1603  }
1604 
1605  switch ( reason )/*lint --e{788}*/
1606  {
1607  case DSDP_CONVERGED:
1608  return 0;
1609 
1610  case DSDP_INFEASIBLE_START:
1611  return 1;
1612 
1613  case DSDP_SMALL_STEPS:
1614  return 2;
1615 
1616  case DSDP_INDEFINITE_SCHUR_MATRIX:
1617  return 2;
1618 
1619  case DSDP_MAX_IT:
1620  return 4;
1621 
1622  case DSDP_NUMERICAL_ERROR:
1623  return 2;
1624 
1625  case DSDP_UPPERBOUND:
1626  return 3;
1627 
1628  case DSDP_USER_TERMINATION:
1629  return 6;
1630 
1631  default:
1632  return 7;
1633  }
1634 }
1635 
1638  SCIP_SDPISOLVER* sdpisolver
1639  )
1640 {
1641  assert( sdpisolver != NULL );
1642  return (SCIPsdpiSolverIsConverged(sdpisolver) && SCIPsdpiSolverIsPrimalFeasible(sdpisolver) && SCIPsdpiSolverIsDualFeasible(sdpisolver));
1643 }
1644 
1648  SCIP_SDPISOLVER* sdpisolver
1649  )
1650 {
1651  assert( sdpisolver != NULL );
1652 
1653  if ( SCIPsdpiSolverIsConverged(sdpisolver) )
1654  return TRUE;
1655 
1656  return FALSE;
1657 }
1658 
1661  SCIP_SDPISOLVER* sdpisolver,
1662  SCIP_Bool* success
1663  )
1664 {/*lint --e{715}*/
1665  SCIPdebugMessage("Not implemented yet\n");
1666  return SCIP_LPERROR;
1667 }
1668 
1671  SCIP_SDPISOLVER* sdpisolver,
1672  SCIP_Real* objval
1673  )
1674 {
1675  assert( sdpisolver != NULL );
1676  assert( objval != NULL );
1677  CHECK_IF_SOLVED( sdpisolver );
1678 
1679  DSDP_CALL( DSDPGetDObjective(sdpisolver->dsdp, objval) );
1680  *objval = -1*(*objval); /*DSDP maximizes instead of minimizing, so the objective values were multiplied by -1 when inserted */
1681 
1682  /* as we didn't add the fixed (lb = ub) variables to dsdp, we have to add their contributions to the objective by hand */
1683  *objval += sdpisolver->fixedvarsobjcontr;
1684 
1685  return SCIP_OKAY;
1686 }
1687 
1693  SCIP_SDPISOLVER* sdpisolver,
1694  SCIP_Real* objval,
1695  SCIP_Real* dualsol,
1696  int* dualsollength
1698  )
1699 {
1700  double* dsdpsol;
1701  int v;
1702  int dsdpnvars;
1703 
1704  assert( sdpisolver != NULL );
1705  assert( dualsollength != NULL );
1706  CHECK_IF_SOLVED( sdpisolver );
1707 
1708  dsdpnvars = sdpisolver->penaltyworbound ? sdpisolver->nactivevars + 1 : sdpisolver->nactivevars; /* in the first case we added r as an explicit var */
1709 
1710  if ( objval != NULL )
1711  {
1712  DSDP_CALL( DSDPGetDObjective(sdpisolver->dsdp, objval) );
1713  *objval *= -1; /* DSDP maximizes instead of minimizing, so the objective values were multiplied by -1 when inserted */
1714 
1715  /* as we didn't add the fixed (lb = ub) variables to dsdp, we have to add their contributions to the objective by hand */
1716  *objval += sdpisolver->fixedvarsobjcontr;
1717  }
1718 
1719  if ( *dualsollength > 0 )
1720  {
1721  assert( dualsol != NULL );
1722  if ( *dualsollength < sdpisolver->nvars )
1723  {
1724  SCIPdebugMessage("The given array in SCIPsdpiSolverGetSol only had length %d, but %d was needed", *dualsollength, sdpisolver->nvars);
1725  *dualsollength = sdpisolver->nvars;
1726 
1727  return SCIP_OKAY;
1728  }
1729 
1730  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &dsdpsol, dsdpnvars) );
1731  DSDP_CALL( DSDPGetY(sdpisolver->dsdp, dsdpsol, dsdpnvars) ); /* last entry needs to be the number of variables, will return an error otherwise */
1732 
1733  /* insert the entries into dualsol, for non-fixed vars we copy those from dsdp, the rest are the saved entries from inserting (they equal lb=ub) */
1734  for (v = 0; v < sdpisolver->nvars; v++)
1735  {
1736  if (sdpisolver->inputtodsdpmapper[v] > -1)
1737  {
1738  /* minus one because the inputtodsdpmapper gives the dsdp indices which start at one, but the array starts at 0 */
1739  dualsol[v] = dsdpsol[sdpisolver->inputtodsdpmapper[v] - 1];
1740  }
1741  else
1742  {
1743  /* this is the value that was saved when inserting, as this variable has lb=ub */
1744  dualsol[v] = sdpisolver->fixedvarsval[(-1 * sdpisolver->inputtodsdpmapper[v]) - 1]; /*lint !e679*/
1745  }
1746  }
1747  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &dsdpsol, dsdpnvars);
1748  }
1749  return SCIP_OKAY;
1750 }
1751 
1760  SCIP_SDPISOLVER* sdpisolver,
1761  SCIP_Real* lbvars,
1762  SCIP_Real* ubvars,
1763  int* arraylength
1765  )
1766 {
1767  double* lbvarsdsdp;
1768  double* ubvarsdsdp;
1769  int i;
1770 
1771  assert( sdpisolver != NULL );
1772  assert( lbvars != NULL );
1773  assert( ubvars != NULL );
1774  assert( arraylength != NULL );
1775  assert( *arraylength >= 0 );
1776  CHECK_IF_SOLVED( sdpisolver );
1777 
1778  /* check if the arrays are long enough */
1779  if ( *arraylength < sdpisolver->nvars )
1780  {
1781  *arraylength = sdpisolver->nvars;
1782  SCIPdebugMessage("Insufficient length of array in SCIPsdpiSolverGetPrimalBoundVars (gave %d, needed %d)\n", *arraylength, sdpisolver->nvars);
1783  return SCIP_OKAY;
1784  }
1785 
1786  /* allocate memory for the arrays given to DSDP */
1787  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &lbvarsdsdp, sdpisolver->nactivevars) );
1788  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &ubvarsdsdp, sdpisolver->nactivevars) );
1789 
1790  /* get the values for the active variables from DSDP */
1791  DSDP_CALL( BConeCopyX(sdpisolver->bcone, lbvarsdsdp, ubvarsdsdp, sdpisolver->nactivevars) );
1792 
1793  /* copy them to the right spots of lbvars & ubvars */
1794  for (i = 0; i < sdpisolver->nvars; i++)
1795  {
1796  if ( sdpisolver->inputtodsdpmapper[i] < 0 )
1797  {
1798  /* if the variable was fixed, it didn't exist in the relaxation, so we set the value to 0
1799  * (as DSDP already uses this value for unbounded vars) */
1800  lbvars[i] = 0;
1801  ubvars[i] = 0;
1802  }
1803  else
1804  {
1805  lbvars[i] = lbvarsdsdp[sdpisolver->inputtodsdpmapper[i] - 1];
1806  ubvars[i] = ubvarsdsdp[sdpisolver->inputtodsdpmapper[i] - 1];
1807  }
1808  }
1809 
1810  /* free allocated memory */
1811  BMSfreeBlockMemoryArrayNull(sdpisolver->blkmem, &ubvarsdsdp, sdpisolver->nactivevars);
1812  BMSfreeBlockMemoryArrayNull(sdpisolver->blkmem, &lbvarsdsdp, sdpisolver->nactivevars);
1813 
1814  return SCIP_OKAY;
1815 }
1816 
1819  SCIP_SDPISOLVER* sdpisolver,
1820  int* iterations
1821  )
1822 {
1823  assert( sdpisolver != NULL );
1824  assert( iterations != NULL );
1825  CHECK_IF_SOLVED( sdpisolver );
1826 
1827  DSDP_CALL( DSDPGetIts(sdpisolver->dsdp, iterations) );
1828  return SCIP_OKAY;
1829 }
1830 
1836 /*
1837  * Numerical Methods
1838  */
1839 
1845  SCIP_SDPISOLVER* sdpisolver
1846  )
1847 {
1848  return 1E+20; /* default infinity from SCIP */
1849 }
1850 
1853  SCIP_SDPISOLVER* sdpisolver,
1854  SCIP_Real val
1855  )
1856 {
1857  return ((val <= -SCIPsdpiSolverInfinity(sdpisolver)) || (val >= SCIPsdpiSolverInfinity(sdpisolver)));
1858 }
1859 
1862  SCIP_SDPISOLVER* sdpisolver
1863  )
1864 {/*lint --e{715}*/
1865  return 1E+10; /* DSDP will start with penalty param 10^10 if called normally */
1866 }
1867 
1870  SCIP_SDPISOLVER* sdpisolver,
1871  SCIP_Real val
1872  )
1873 {
1874  return ((val <= -SCIPsdpiSolverMaxPenParam(sdpisolver)) || (val >= SCIPsdpiSolverMaxPenParam(sdpisolver)));
1875 }
1876 
1879  SCIP_SDPISOLVER* sdpisolver,
1880  SCIP_SDPPARAM type,
1881  SCIP_Real* dval
1882  )
1883 {
1884  assert( sdpisolver != NULL );
1885  assert( dval != NULL );
1886 
1887  switch( type )/*lint --e{788}*/
1888  {
1889  case SCIP_SDPPAR_EPSILON:
1890  *dval = sdpisolver->epsilon;
1891  break;
1892  case SCIP_SDPPAR_FEASTOL:
1893  *dval = sdpisolver->feastol;
1894  break;
1895  case SCIP_SDPPAR_OBJLIMIT:
1896  *dval = sdpisolver->objlimit;
1897  break;
1898  default:
1899  return SCIP_PARAMETERUNKNOWN;
1900  }
1901 
1902  return SCIP_OKAY;
1903 }
1904 
1907  SCIP_SDPISOLVER* sdpisolver,
1908  SCIP_SDPPARAM type,
1909  SCIP_Real dval
1910  )
1911 {
1912  assert( sdpisolver != NULL );
1913 
1914  switch( type )/*lint --e{788}*/
1915  {
1916  case SCIP_SDPPAR_EPSILON:
1917  sdpisolver->epsilon = dval;
1918  SCIPdebugMessage("Setting sdpisolver epsilon to %f.\n", dval);
1919  break;
1920  case SCIP_SDPPAR_FEASTOL:
1921  sdpisolver->feastol = dval;
1922  SCIPdebugMessage("Setting sdpisolver feastol to %f.\n", dval);
1923  break;
1924  case SCIP_SDPPAR_OBJLIMIT:
1925  /* DSDP only allows to set a dual bound, but as we want to solve the dual problem in DSDP, we would need to set a primal bound, which doesn't exist in
1926  * DSDP, so we can't do anything in this case. */
1927  SCIPdebugMessage("Setting sdpisolver objlimit to %f.\n", dval);
1928  sdpisolver->objlimit = dval;
1929  break;
1930  default:
1931  return SCIP_PARAMETERUNKNOWN;
1932  }
1933 
1934  return SCIP_OKAY;
1935 }
1936 
1939  SCIP_SDPISOLVER* sdpisolver,
1940  SCIP_SDPPARAM type,
1941  int* ival
1942  )
1943 {
1944  assert( sdpisolver != NULL );
1945 
1946  switch( type )/*lint --e{788}*/
1947  {
1948  case SCIP_SDPPAR_SDPINFO:
1949  *ival = sdpisolver->sdpinfo;
1950  SCIPdebugMessage("Getting sdpisolver information output (%d).\n", *ival);
1951  break;
1952  default:
1953  return SCIP_PARAMETERUNKNOWN;
1954  }
1955 
1956  return SCIP_OKAY;
1957 }
1958 
1961  SCIP_SDPISOLVER* sdpisolver,
1962  SCIP_SDPPARAM type,
1963  int ival
1964  )
1965 {
1966  assert( sdpisolver != NULL );
1967 
1968  switch( type )/*lint --e{788}*/
1969  {
1970  case SCIP_SDPPAR_SDPINFO:
1971  sdpisolver->sdpinfo = (SCIP_Bool) ival;
1972  SCIPdebugMessage("Setting sdpisolver information output (%d).\n", ival);
1973  break;
1974  default:
1975  return SCIP_PARAMETERUNKNOWN;
1976  }
1977 
1978  return SCIP_OKAY;
1979 }
1980 
1986 /*
1987  * File Interface Methods
1988  */
1989 
1995  SCIP_SDPISOLVER* sdpisolver,
1996  const char* fname
1997  )
1998 {/*lint --e{715}*/
1999  SCIPdebugMessage("Not implemented yet\n");
2000  return SCIP_LPERROR;
2001 }
2002 
2005  SCIP_SDPISOLVER* sdpisolver,
2006  const char* fname
2007  )
2008 {/*lint --e{715}*/
2009  SCIPdebugMessage("Not implemented yet\n");
2010  return SCIP_LPERROR;
2011 }
2012 
SCIP_Bool SCIPsdpiSolverIsConverged(SCIP_SDPISOLVER *sdpisolver)
int SCIPsdpiSolverGetInternalStatus(SCIP_SDPISOLVER *sdpisolver)
SCIP_Bool SCIPsdpiSolverIsPrimalFeasible(SCIP_SDPISOLVER *sdpisolver)
SCIP_Bool SCIPsdpiSolverIsPrimalUnbounded(SCIP_SDPISOLVER *sdpisolver)
SCIP_Real SCIPsdpiSolverMaxPenParam(SCIP_SDPISOLVER *sdpisolver)
SCIP_Bool SCIPsdpiSolverIsObjlimExc(SCIP_SDPISOLVER *sdpisolver)
SCIP_Bool SCIPsdpiSolverIsDualFeasible(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverGetSol(SCIP_SDPISOLVER *sdpisolver, SCIP_Real *objval, SCIP_Real *dualsol, int *dualsollength)
SCIP_Real SCIPsdpiSolverInfinity(SCIP_SDPISOLVER *sdpisolver)
static SCIP_Bool isFixed(SCIP_SDPISOLVER *sdpisolver, SCIP_Real lb, SCIP_Real ub)
const char * SCIPsdpiSolverGetSolverDesc(void)
#define CHECK_IF_SOLVED_BOOL(sdpisolver)
SCIP_RETCODE SCIPsdpiSolverGetIntpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, int *ival)
#define DSDP_CALL_BOOL(x)
interface methods for specific SDP solvers
SCIP_RETCODE SCIPsdpiSolverResetCounter(SCIP_SDPISOLVER *sdpisolver)
void * SCIPsdpiSolverGetSolverPointer(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverGetRealpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, SCIP_Real *dval)
SCIP_Bool SCIPsdpiSolverIsInfinity(SCIP_SDPISOLVER *sdpisolver, SCIP_Real val)
SCIP_RETCODE SCIPsdpiSolverSetRealpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, SCIP_Real dval)
SCIP_Bool SCIPsdpiSolverIsPrimalInfeasible(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverGetSolFeasibility(SCIP_SDPISOLVER *sdpisolver, SCIP_Bool *primalfeasible, SCIP_Bool *dualfeasible)
SCIP_RETCODE SCIPsdpiSolverFree(SCIP_SDPISOLVER **sdpisolver)
#define DSDP_CALLM(x)
#define CHECK_IF_SOLVED(sdpisolver)
static void sortColRow(int *row, int *col, SCIP_Real *val, int length)
SCIP_Bool SCIPsdpiSolverIsDualUnbounded(SCIP_SDPISOLVER *sdpisolver)
SCIP_Bool SCIPsdpiSolverFeasibilityKnown(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverLoadAndSolve(SCIP_SDPISOLVER *sdpisolver, int nvars, SCIP_Real *obj, SCIP_Real *lb, SCIP_Real *ub, int nsdpblocks, int *sdpblocksizes, int *sdpnblockvars, int sdpconstnnonz, int *sdpconstnblocknonz, int **sdpconstrow, int **sdpconstcol, SCIP_Real **sdpconstval, int sdpnnonz, int **sdpnblockvarnonz, int **sdpvar, int ***sdprow, int ***sdpcol, SCIP_Real ***sdpval, int **indchanges, int *nremovedinds, int *blockindchanges, int nremovedblocks, int nlpcons, int noldlpcons, SCIP_Real *lplhs, SCIP_Real *lprhs, int *rownactivevars, int lpnnonz, int *lprow, int *lpcol, SCIP_Real *lpval, SCIP_Real *start)
SCIP_RETCODE SCIPsdpiSolverGetObjval(SCIP_SDPISOLVER *sdpisolver, SCIP_Real *objval)
#define BMS_CALL(x)
const char * SCIPsdpiSolverGetSolverName(void)
SCIP_Bool SCIPsdpiSolverIsGEMaxPenParam(SCIP_SDPISOLVER *sdpisolver, SCIP_Real val)
SCIP_RETCODE SCIPsdpiSolverSetIntpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, int ival)
SCIP_Bool SCIPsdpiSolverIsIterlimExc(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverGetPrimalBoundVars(SCIP_SDPISOLVER *sdpisolver, SCIP_Real *lbvars, SCIP_Real *ubvars, int *arraylength)
SCIP_RETCODE SCIPsdpiSolverIncreaseCounter(SCIP_SDPISOLVER *sdpisolver)
static int compLowerTriangPos(int i, int j)
SCIP_Bool SCIPsdpiSolverIsTimelimExc(SCIP_SDPISOLVER *sdpisolver)
SCIP_Bool SCIPsdpiSolverWasSolved(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverCreate(SCIP_SDPISOLVER **sdpisolver, SCIP_MESSAGEHDLR *messagehdlr, BMS_BLKMEM *blkmem)
SCIP_RETCODE SCIPsdpiSolverIgnoreInstability(SCIP_SDPISOLVER *sdpisolver, SCIP_Bool *success)
struct SCIP_SDPiSolver SCIP_SDPISOLVER
Definition: sdpisolver.h:71
SCIP_Bool SCIPsdpiSolverIsOptimal(SCIP_SDPISOLVER *sdpisolver)
SCIP_Bool SCIPsdpiSolverIsAcceptable(SCIP_SDPISOLVER *sdpisolver)
enum SCIP_SDPParam SCIP_SDPPARAM
Definition: type_sdpi.h:63
SCIP_Bool SCIPsdpiSolverIsDualInfeasible(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverLoadAndSolveWithPenalty(SCIP_SDPISOLVER *sdpisolver, SCIP_Real penaltyparam, SCIP_Bool withobj, SCIP_Bool rbound, int nvars, SCIP_Real *obj, SCIP_Real *lb, SCIP_Real *ub, int nsdpblocks, int *sdpblocksizes, int *sdpnblockvars, int sdpconstnnonz, int *sdpconstnblocknonz, int **sdpconstrow, int **sdpconstcol, SCIP_Real **sdpconstval, int sdpnnonz, int **sdpnblockvarnonz, int **sdpvar, int ***sdprow, int ***sdpcol, SCIP_Real ***sdpval, int **indchanges, int *nremovedinds, int *blockindchanges, int nremovedblocks, int nlpcons, int noldlpcons, SCIP_Real *lplhs, SCIP_Real *lprhs, int *rownactivevars, int lpnnonz, int *lprow, int *lpcol, SCIP_Real *lpval, SCIP_Real *start, SCIP_Bool *feasorig)
SCIP_RETCODE SCIPsdpiSolverWriteSDP(SCIP_SDPISOLVER *sdpisolver, const char *fname)
SCIP_RETCODE SCIPsdpiSolverReadSDP(SCIP_SDPISOLVER *sdpisolver, const char *fname)
SCIP_RETCODE SCIPsdpiSolverGetIterations(SCIP_SDPISOLVER *sdpisolver, int *iterations)
#define DSDP_CALL(x)