SCIP-SDP  4.0.0
sdpsolchecker.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 programs based on SCIP. */
5 /* */
6 /* Copyright (C) 2011-2013 Discrete Optimization, TU Darmstadt */
7 /* EDOM, FAU Erlangen-Nürnberg */
8 /* 2014-2021 Discrete Optimization, TU Darmstadt */
9 /* */
10 /* */
11 /* This program is free software; you can redistribute it and/or */
12 /* modify it under the terms of the GNU Lesser General Public License */
13 /* as published by the Free Software Foundation; either version 3 */
14 /* of the License, or (at your option) any later version. */
15 /* */
16 /* This program is distributed in the hope that it will be useful, */
17 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
18 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
19 /* GNU Lesser General Public License for more details. */
20 /* */
21 /* You should have received a copy of the GNU Lesser General Public License */
22 /* along with this program; if not, write to the Free Software */
23 /* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.*/
24 /* */
25 /* */
26 /* Based on SCIP - Solving Constraint Integer Programs */
27 /* Copyright (C) 2002-2021 Zuse Institute Berlin */
28 /* SCIP is distributed under the terms of the SCIP Academic Licence, */
29 /* see file COPYING in the SCIP distribution. */
30 /* */
31 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
32 
39 /* #define SCIP_DEBUG */
40 /* #define SCIP_MORE_DEBUG /\* shows all values of variables and constraint and their violations *\/ */
41 
42 #include "sdpi/sdpsolchecker.h"
43 #include "sdpi/lapack_interface.h"
44 #include "scip/pub_message.h" /* for debug and error message */
45 
47 #define BMS_CALL(x) do \
48  { \
49  if( NULL == (x) ) \
50  { \
51  SCIPerrorMessage("No memory in function call.\n"); \
52  return SCIP_NOMEMORY; \
53  } \
54  } \
55  while( FALSE )
56 
57 #define INF 1e+16
58 
62 SCIP_RETCODE SCIPsdpSolcheckerCheck(
63  BMS_BUFMEM* bufmem,
64  int nvars,
65  SCIP_Real* lb,
66  SCIP_Real* ub,
67  int nsdpblocks,
68  int* sdpblocksizes,
69  int* sdpnblockvars,
70  int sdpconstnnonz,
71  int* sdpconstnblocknonz,
73  int** sdpconstrow,
74  int** sdpconstcol,
75  SCIP_Real** sdpconstval,
76  int sdpnnonz,
77  int** sdpnblockvarnonz,
79  int** sdpvar,
81  int*** sdprow,
82  int*** sdpcol,
83  SCIP_Real*** sdpval,
84  int** indchanges,
86  int* nremovedinds,
87  int* blockindchanges,
88  int nlpcons,
89  SCIP_Real* lplhs,
90  SCIP_Real* lprhs,
91  int lpnnonz,
92  int* lprow,
93  int* lpcol,
94  SCIP_Real* lpval,
95  SCIP_Real* solvector,
96  SCIP_Real feastol,
97  SCIP_Real epsilon,
98  SCIP_Bool* infeasible
99  )
100 {/*lint --e{818}*/
101  int i;
102  int j;
103  int b;
104  int v;
105  int ind;
106 
107  assert( bufmem != NULL );
108  assert( lb != NULL );
109  assert( ub != NULL );
110  assert( nsdpblocks >= 0 );
111  assert( nsdpblocks == 0 || sdpblocksizes != NULL );
112  assert( nsdpblocks == 0 || sdpnblockvars != NULL );
113  assert( sdpconstnnonz >= 0 );
114  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstnblocknonz != NULL );
115  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstrow != NULL );
116  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstcol != NULL );
117  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstval != NULL );
118  assert( sdpnnonz >= 0 );
119  assert( nsdpblocks == 0 || sdpnblockvarnonz != NULL );
120  assert( nsdpblocks == 0 || sdpvar != NULL );
121  assert( nsdpblocks == 0 || sdprow != NULL );
122  assert( nsdpblocks == 0 || sdpcol != NULL );
123  assert( nsdpblocks == 0 || sdpval != NULL );
124  assert( nsdpblocks == 0 || indchanges != NULL );
125  assert( nsdpblocks == 0 || nremovedinds != NULL );
126  assert( nsdpblocks == 0 || blockindchanges != NULL );
127  assert( nlpcons >= 0 );
128  assert( nlpcons == 0 || lplhs != NULL );
129  assert( nlpcons == 0 || lprhs != NULL );
130  assert( lpnnonz >= 0 );
131  assert( nlpcons == 0 || lprow != NULL );
132  assert( nlpcons == 0 || lpcol != NULL );
133  assert( nlpcons == 0 || lpval != NULL );
134  assert( solvector != NULL );
135  assert( feastol >= 0 );
136  assert( infeasible != NULL );
137 
138  /* check variable bounds */
139  for (i = 0; i < nvars; i++)
140  {
141  if ( solvector[i] < lb[i] - feastol || solvector[i] > ub[i] + feastol )
142  {
143  SCIPdebugMessage("solution found infeasible (feastol=%g) for dual variable bounds: x[%d] = %g <|= [%g, %g]\n",
144  feastol, i, solvector[i], lb[i], ub[i]);
145  *infeasible = TRUE;
146  return SCIP_OKAY;
147  }
148  }
149 
150  /* check linear constraints (since DSDP sorts the lp-arrays by cols, we cannot expect them to be sorted) */
151  if ( nlpcons > 0 )
152  {
153  SCIP_Real* lpconsvals;
154 
155  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &lpconsvals, nlpcons) );
156 
157  /* initialize all rows with zero */
158  for (i = 0; i < nlpcons; i++)
159  lpconsvals[i] = 0;
160 
161  /* compute the values of all rows */
162  for (i = 0; i < lpnnonz; i++)
163  {
164  if ( lb[lpcol[i]] < ub[lpcol[i]] - epsilon ) /* fixed variables are already included in lhs/rhs */
165  lpconsvals[lprow[i]] += solvector[lpcol[i]] * lpval[i];
166  }
167 
168  /* check all active constraints for feasibility */
169  ind = 0; /* used to iterate over active constraints */
170  for (i = 0; i < nlpcons; i++)
171  {
172  if ( lpconsvals[i] < lplhs[ind] - feastol || lpconsvals[i] > lprhs[ind] + feastol)
173  {
174  SCIPdebugMessage("solution found infeasible (feastol=%g) for dual lp constraint: LP-%d = %g <|= [%g,%g]\n",
175  feastol, i, lpconsvals[i], lplhs[ind], lprhs[ind]);
176  BMSfreeBufferMemoryArray(bufmem, &lpconsvals);
177  *infeasible = TRUE;
178  return SCIP_OKAY;
179  }
180 
181  ind++;
182  }
183  BMSfreeBufferMemoryArray(bufmem, &lpconsvals);
184  }
185 
186  /* check sdp constraints */
187  if ( nsdpblocks > 0 )
188  {
189  SCIP_Real* fullsdpmatrix;
190  SCIP_Real eigenvalue;
191  int maxblocksize = 0;
192 
193  /* allocate memory */
194  if ( nsdpblocks == 1 )
195  maxblocksize = sdpblocksizes[0] - nremovedinds[0];
196  else
197  {
198  /* calculate maximum size of any SDP block to not have to reallocate memory in between */
199  for (b = 0; b < nsdpblocks; b++)
200  {
201  if ( (sdpblocksizes[b] - nremovedinds[b]) > maxblocksize )
202  maxblocksize = sdpblocksizes[b] - nremovedinds[b];
203  }
204  }
205 
206  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &fullsdpmatrix, maxblocksize * maxblocksize) );/*lint !e647*/
207 
208  for (b = 0; b < nsdpblocks; b++)
209  {
210  if ( blockindchanges[b] > -1 )
211  {
212  /* initialize lower triangular part of fullsdpmatrix with zero */
213  for (i = 0; i < sdpblocksizes[b] - nremovedinds[b]; i++)
214  {
215  for (j = 0; j <= i; j++)
216  {
217  fullsdpmatrix[i * (sdpblocksizes[b] - nremovedinds[b]) + j] = 0.0; /*lint !e679*/
218  }
219  }
220 
221  /* iterate over all non-fixed variables and add the corresponding nonzeros */
222  for (v = 0; v < sdpnblockvars[b]; v++)
223  {
224  if ( lb[sdpvar[b][v]] < ub[sdpvar[b][v]] - epsilon )
225  {
226  for (i = 0; i < sdpnblockvarnonz[b][v]; i++)
227  {
228  fullsdpmatrix[((sdprow[b][v][i] - indchanges[b][sdprow[b][v][i]]) * (sdpblocksizes[b] - nremovedinds[b])) +
229  sdpcol[b][v][i] - indchanges[b][sdpcol[b][v][i]]] += solvector[sdpvar[b][v]] * sdpval[b][v][i];/*lint !e679*/
230  }
231  }
232  }
233 
234  /* add constant matrix */
235  if ( sdpconstnblocknonz != NULL )
236  {
237  for (i = 0; i < sdpconstnblocknonz[b]; i++)
238  {
239  fullsdpmatrix[((sdpconstrow[b][i] - indchanges[b][sdpconstrow[b][i]]) * (sdpblocksizes[b] - nremovedinds[b])) +
240  sdpconstcol[b][i] - indchanges[b][sdpconstcol[b][i]]] -= sdpconstval[b][i];/*lint !e679*/
241  }
242  }
243 
244  /* extend to full symmetric matrix for LAPACK */
245  for (i = 0; i < sdpblocksizes[b] - nremovedinds[b]; i++)
246  {
247  for (j = 0; j < i; j++)
248  {
249  fullsdpmatrix[j * (sdpblocksizes[b] - nremovedinds[b]) + i] = fullsdpmatrix[i * (sdpblocksizes[b] - nremovedinds[b]) + j];/*lint !e679*/
250  }
251  }
252 
253  /* compute smallest eigenvalue using LAPACK */
254  SCIP_CALL( SCIPlapackComputeIthEigenvalue(bufmem, FALSE, sdpblocksizes[b] - nremovedinds[b], fullsdpmatrix, 1, &eigenvalue, NULL) );
255 
256  if ( eigenvalue < - feastol )
257  {
258  SCIPdebugMessage("solution found infeasible (feastol=%g) for dual sdp constraint %d, smallest eigenvector %.10g\n",
259  feastol, b, eigenvalue);
260  BMSfreeBufferMemoryArray(bufmem, &fullsdpmatrix);
261  *infeasible = TRUE;
262  return SCIP_OKAY;
263  }
264  }
265  }
266 
267  BMSfreeBufferMemoryArray(bufmem, &fullsdpmatrix);
268  }
269 
270  *infeasible = FALSE;
271  return SCIP_OKAY;
272 }
273 
274 
281 SCIP_RETCODE SCIPsdpSolcheckerCheckAndGetViolDual(
282  BMS_BUFMEM* bufmem,
283  int nvars,
284  SCIP_Real* lb,
285  SCIP_Real* ub,
286  int nsdpblocks,
287  int* sdpblocksizes,
288  int* sdpnblockvars,
289  int sdpconstnnonz,
290  int* sdpconstnblocknonz,
292  int** sdpconstrow,
293  int** sdpconstcol,
294  SCIP_Real** sdpconstval,
295  int sdpnnonz,
296  int** sdpnblockvarnonz,
298  int** sdpvar,
300  int*** sdprow,
301  int*** sdpcol,
302  SCIP_Real*** sdpval,
303  int** indchanges,
305  int* nremovedinds,
306  int* blockindchanges,
307  int nlpcons,
308  SCIP_Real* lplhs,
309  SCIP_Real* lprhs,
310  int lpnnonz,
311  int* lprow,
312  int* lpcol,
313  SCIP_Real* lpval,
314  SCIP_Real* solvector,
315  SCIP_Real feastol,
316  SCIP_Real epsilon,
317  SCIP_Real* maxabsviolbnds,
318  SCIP_Real* sumabsviolbnds,
319  SCIP_Real* maxabsviolcons,
320  SCIP_Real* sumabsviolcons,
321  SCIP_Real* maxabsviolsdp,
322  SCIP_Real* sumabsviolsdp,
323  SCIP_Bool* infeasible
324  )
325 {/*lint --e{818}*/
326  int i;
327  int j;
328  int b;
329  int v;
330  int ind;
331  SCIP_Real viol;
332 
333  assert( bufmem != NULL );
334  assert( lb != NULL );
335  assert( ub != NULL );
336  assert( nsdpblocks >= 0 );
337  assert( nsdpblocks == 0 || sdpblocksizes != NULL );
338  assert( nsdpblocks == 0 || sdpnblockvars != NULL );
339  assert( sdpconstnnonz >= 0 );
340  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstnblocknonz != NULL );
341  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstrow != NULL );
342  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstcol != NULL );
343  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstval != NULL );
344  assert( sdpnnonz >= 0 );
345  assert( nsdpblocks == 0 || sdpnblockvarnonz != NULL );
346  assert( nsdpblocks == 0 || sdpvar != NULL );
347  assert( nsdpblocks == 0 || sdprow != NULL );
348  assert( nsdpblocks == 0 || sdpcol != NULL );
349  assert( nsdpblocks == 0 || sdpval != NULL );
350  assert( nsdpblocks == 0 || indchanges != NULL );
351  assert( nsdpblocks == 0 || nremovedinds != NULL );
352  assert( nsdpblocks == 0 || blockindchanges != NULL );
353  assert( nlpcons >= 0 );
354  assert( nlpcons == 0 || lplhs != NULL );
355  assert( nlpcons == 0 || lprhs != NULL );
356  assert( lpnnonz >= 0 );
357  assert( nlpcons == 0 || lprow != NULL );
358  assert( nlpcons == 0 || lpcol != NULL );
359  assert( nlpcons == 0 || lpval != NULL );
360  assert( solvector != NULL );
361  assert( feastol >= 0 );
362  assert( maxabsviolbnds != NULL );
363  assert( sumabsviolbnds != NULL );
364  assert( maxabsviolcons != NULL );
365  assert( sumabsviolcons != NULL );
366  assert( maxabsviolsdp != NULL );
367  assert( sumabsviolsdp != NULL );
368  assert( infeasible != NULL );
369 
370  /* initialize violations with 0.0 */
371  *maxabsviolbnds = 0.0;
372  *sumabsviolbnds = 0.0;
373  *maxabsviolcons = 0.0;
374  *sumabsviolcons = 0.0;
375  *maxabsviolsdp = 0.0;
376  *sumabsviolsdp = 0.0;
377 
378  *infeasible = FALSE;
379 
380  /* check variable bounds */
381  for (i = 0; i < nvars; i++)
382  {
383  viol = MAX3(lb[i] - solvector[i], solvector[i] - ub[i], 0.0);
384  *sumabsviolbnds += viol;
385 
386 #ifdef SCIP_MORE_DEBUG
387  SCIPdebugMessage("Dual variable %d: lb = %g, ub = %g, val = %g, viol = %g\n", i, lb[i], ub[i], solvector[i], viol);
388 #endif
389 
390  if ( viol > *maxabsviolbnds )
391  *maxabsviolbnds = viol;
392 
393  if ( solvector[i] < lb[i] - feastol || solvector[i] > ub[i] + feastol )
394  {
395  SCIPdebugMessage("solution found infeasible (feastol=%f) for dual variable bounds: x[%d] = %g <|= [%g, %g]\n",
396  feastol, i, solvector[i], lb[i], ub[i]);
397  *infeasible = TRUE;
398  }
399  }
400 
401  /* check linear constraints (since DSDP sorts the lp-arrays by cols, we cannot expect them to be sorted) */
402  if ( nlpcons > 0 )
403  {
404  SCIP_Real* lpconsvals;
405 
406  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &lpconsvals, nlpcons) );
407 
408  /* initialize all rows with zero */
409  for (i = 0; i < nlpcons; i++)
410  lpconsvals[i] = 0;
411 
412  /* compute the values of all rows */
413  for (i = 0; i < lpnnonz; i++)
414  {
415  if ( lb[lpcol[i]] < ub[lpcol[i]] - epsilon ) /* fixed variables are already included in lhs/rhs */
416  lpconsvals[lprow[i]] += solvector[lpcol[i]] * lpval[i];
417  }
418 
419  /* check all active constraints for feasibility */
420  ind = 0; /* used to iterate over active constraints */
421  for (i = 0; i < nlpcons; i++)
422  {
423  viol = MAX3(lplhs[ind] - lpconsvals[i], lpconsvals[i] - lprhs[ind], 0.0);
424  *sumabsviolcons += viol;
425 
426 #ifdef SCIP_MORE_DEBUG
427  SCIPdebugMessage("Dual linear constraint %d: lhs = %g, rhs = %g, val = %.15g, viol = %.15g\n", i, lplhs[i], lprhs[i], lpconsvals[i], viol);
428 #endif
429 
430  if ( viol > *maxabsviolcons )
431  *maxabsviolcons = viol;
432 
433  if ( lpconsvals[i] < lplhs[ind] - feastol || lpconsvals[i] > lprhs[ind] + feastol)
434  {
435  SCIPdebugMessage("solution found infeasible (feastol=%g) for dual lp constraint: LP-%d = %g <|= [%g,%g].\n",
436  feastol, i, lpconsvals[i], lplhs[ind], lprhs[ind]);
437  BMSfreeBufferMemoryArray(bufmem, &lpconsvals);
438  *infeasible = TRUE;
439  }
440 
441  ind++;
442  }
443  BMSfreeBufferMemoryArray(bufmem, &lpconsvals);
444  }
445 
446  /* check sdp constraints */
447  if ( nsdpblocks > 0 )
448  {
449  SCIP_Real* fullsdpmatrix;
450  SCIP_Real eigenvalue;
451  int maxblocksize = 0;
452 
453  /* allocate memory */
454  if ( nsdpblocks == 1 )
455  maxblocksize = sdpblocksizes[0] - nremovedinds[0];
456  else
457  {
458  /* calculate maximum size of any SDP block to not have to reallocate memory in between */
459  for (b = 0; b < nsdpblocks; b++)
460  {
461  if ( (sdpblocksizes[b] - nremovedinds[b]) > maxblocksize )
462  maxblocksize = sdpblocksizes[b] - nremovedinds[b];
463  }
464  }
465 
466  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &fullsdpmatrix, maxblocksize * maxblocksize) );/*lint !e647*/
467 
468  for (b = 0; b < nsdpblocks; b++)
469  {
470  if ( blockindchanges[b] > -1 )
471  {
472  /* initialize lower triangular part of fullsdpmatrix with zero */
473  for (i = 0; i < sdpblocksizes[b] - nremovedinds[b]; i++)
474  {
475  for (j = 0; j <= i; j++)
476  {
477  fullsdpmatrix[i * (sdpblocksizes[b] - nremovedinds[b]) + j] = 0.0; /*lint !e679*/
478  }
479  }
480 
481  /* iterate over all non-fixed variables and add the corresponding nonzeros */
482  for (v = 0; v < sdpnblockvars[b]; v++)
483  {
484  if ( lb[sdpvar[b][v]] < ub[sdpvar[b][v]] - epsilon )
485  {
486  for (i = 0; i < sdpnblockvarnonz[b][v]; i++)
487  {
488  fullsdpmatrix[((sdprow[b][v][i] - indchanges[b][sdprow[b][v][i]]) * (sdpblocksizes[b] - nremovedinds[b])) +
489  sdpcol[b][v][i] - indchanges[b][sdpcol[b][v][i]]] += solvector[sdpvar[b][v]] * sdpval[b][v][i];/*lint !e679*/
490  }
491  }
492  }
493 
494  /* add constant matrix */
495  if ( sdpconstnblocknonz != NULL )
496  {
497  for (i = 0; i < sdpconstnblocknonz[b]; i++)
498  {
499  fullsdpmatrix[((sdpconstrow[b][i] - indchanges[b][sdpconstrow[b][i]]) * (sdpblocksizes[b] - nremovedinds[b])) +
500  sdpconstcol[b][i] - indchanges[b][sdpconstcol[b][i]]] -= sdpconstval[b][i];/*lint !e679*/
501  }
502  }
503 
504  /* extend to full symmetric matrix for LAPACK */
505  for (i = 0; i < sdpblocksizes[b] - nremovedinds[b]; i++)
506  {
507  for (j = 0; j < i; j++)
508  {
509  fullsdpmatrix[j * (sdpblocksizes[b] - nremovedinds[b]) + i] = fullsdpmatrix[i * (sdpblocksizes[b] - nremovedinds[b]) + j];/*lint !e679*/
510  }
511  }
512 
513  /* compute smallest eigenvalue using LAPACK */
514  SCIP_CALL( SCIPlapackComputeIthEigenvalue(bufmem, FALSE, sdpblocksizes[b] - nremovedinds[b], fullsdpmatrix, 1, &eigenvalue, NULL) );
515 
516  viol = MAX(-eigenvalue, 0.0);
517  *sumabsviolsdp += viol;
518 
519 #ifdef SCIP_MORE_DEBUG
520  SCIPdebugMessage("Dual SDP constraint %d: lambda_min = %.15g, viol = %.15g\n", b, eigenvalue, viol);
521 #endif
522 
523  if ( viol > *maxabsviolsdp )
524  *maxabsviolsdp = viol;
525 
526  if ( eigenvalue < - feastol )
527  {
528  SCIPdebugMessage("solution found infeasible (feastol=%.10g) for dual SDP constraint %d, smallest eigenvector %.10g\n",
529  feastol, b, eigenvalue);
530  BMSfreeBufferMemoryArray(bufmem, &fullsdpmatrix);
531  *infeasible = TRUE;
532  }
533  }
534  }
535 
536  BMSfreeBufferMemoryArray(bufmem, &fullsdpmatrix);
537  }
538 
539  return SCIP_OKAY;
540 }
541 
542 
550 SCIP_RETCODE SCIPsdpSolcheckerCheckAndGetViolPrimal(
551  BMS_BUFMEM* bufmem,
552  int nvars,
553  SCIP_Real* obj,
554  SCIP_Real* lb,
555  SCIP_Real* ub,
556  int* inputtomosekmapper,
559  int nsdpblocks,
560  int* sdpblocksizes,
561  int* sdpnblockvars,
562  int sdpconstnnonz,
563  int* sdpconstnblocknonz,
565  int** sdpconstrow,
566  int** sdpconstcol,
567  SCIP_Real** sdpconstval,
568  int sdpnnonz,
569  int** sdpnblockvarnonz,
571  int** sdpvar,
573  int*** sdprow,
574  int*** sdpcol,
575  SCIP_Real*** sdpval,
576  int** indchanges,
578  int* nremovedinds,
579  int* blockindchanges,
580  int nremovedblocks,
581  int nlpcons,
582  SCIP_Real* lplhs,
583  SCIP_Real* lprhs,
584  int lpnnonz,
585  int* lprow,
586  int* lpcol,
587  SCIP_Real* lpval,
588  SCIP_Real* solvector,
589  SCIP_Real** solmatrices,
590  SCIP_Real feastol,
591  SCIP_Real epsilon,
592  SCIP_Real* maxabsviolbnds,
593  SCIP_Real* sumabsviolbnds,
594  SCIP_Real* maxabsviolcons,
595  SCIP_Real* sumabsviolcons,
596  SCIP_Real* maxabsviolsdp,
597  SCIP_Real* sumabsviolsdp,
598  SCIP_Bool* infeasible
599  )
600 {/*lint --e{818}*/
601  int i;
602  int j;
603  int b;
604  int v;
605  SCIP_Real viol;
606  int nactivevars;
607  int nvarbounds;
608  int nfixedvars;
609  SCIP_Real* objcoefs;
610  int* varboundpos;
611  int nlpvars;
612  int nprimalvars;
613  int nprimalmatrixvars;
614  int maxblocksize;
615  int* mosekblocksizes;
616  int nprimallpcons;
617  int blocksize;
618  int mosekind;
619  int mosekrow;
620  int mosekcol;
621  int blockvar;
622  int k;
623  int nnonz;
624  int c;
625 
626  assert( bufmem != NULL );
627  assert( lb != NULL );
628  assert( ub != NULL );
629  assert( nsdpblocks >= 0 );
630  assert( nsdpblocks == 0 || sdpblocksizes != NULL );
631  assert( nsdpblocks == 0 || sdpnblockvars != NULL );
632  assert( sdpconstnnonz >= 0 );
633  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstnblocknonz != NULL );
634  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstrow != NULL );
635  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstcol != NULL );
636  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstval != NULL );
637  assert( sdpnnonz >= 0 );
638  assert( nsdpblocks == 0 || sdpnblockvarnonz != NULL );
639  assert( nsdpblocks == 0 || sdpvar != NULL );
640  assert( nsdpblocks == 0 || sdprow != NULL );
641  assert( nsdpblocks == 0 || sdpcol != NULL );
642  assert( nsdpblocks == 0 || sdpval != NULL );
643  assert( nsdpblocks == 0 || indchanges != NULL );
644  assert( nsdpblocks == 0 || nremovedinds != NULL );
645  assert( nsdpblocks == 0 || blockindchanges != NULL );
646  assert( 0 <= nremovedblocks && nremovedblocks <= nsdpblocks );
647  assert( nlpcons >= 0 );
648  assert( nlpcons == 0 || lplhs != NULL );
649  assert( nlpcons == 0 || lprhs != NULL );
650  assert( lpnnonz >= 0 );
651  assert( nlpcons == 0 || lprow != NULL );
652  assert( nlpcons == 0 || lpcol != NULL );
653  assert( nlpcons == 0 || lpval != NULL );
654  assert( solvector != NULL );
655  assert( feastol >= 0 );
656  assert( maxabsviolbnds != NULL );
657  assert( sumabsviolbnds != NULL );
658  assert( maxabsviolcons != NULL );
659  assert( sumabsviolcons != NULL );
660  assert( maxabsviolsdp != NULL );
661  assert( sumabsviolsdp != NULL );
662  assert( infeasible != NULL );
663 
664  /* initialize violations with 0.0 */
665  *maxabsviolbnds = 0.0;
666  *sumabsviolbnds = 0.0;
667  *maxabsviolcons = 0.0;
668  *sumabsviolcons = 0.0;
669  *maxabsviolsdp = 0.0;
670  *sumabsviolsdp = 0.0;
671 
672  *infeasible = FALSE;
673 
674  nactivevars = 0;
675  nvarbounds = 0;
676  nfixedvars = 0;
677 
678  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &objcoefs, nvars) );
679  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &varboundpos, 2 * nvars) );
680 
681  /* find fixed variables in dual problem */
682  for (i = 0; i < nvars; i++)
683  {
684  if ( ub[i] - lb[i] <= epsilon )
685  nfixedvars++;
686  else
687  {
688  objcoefs[nactivevars] = obj[i];
689 
690  if ( lb[i] > -INF )
691  varboundpos[nvarbounds++] = -(nactivevars + 1); /* negative sign means lower bound */
692 
693  if ( ub[i] < INF )
694  varboundpos[nvarbounds++] = +(nactivevars + 1); /* positive sign means upper bound */
695 
696  nactivevars++;
697  }
698  }
699  assert( nactivevars + nfixedvars == nvars );
700 
701  /* determine total number of sides in LP-constraints */
702  nlpvars = 0;
703  if ( nlpcons > 0 )
704  {
705  for (i = 0; i < nlpcons; i++)
706  {
707  if ( lplhs[i] > - INF )
708  ++nlpvars;
709 
710  if ( lprhs[i] < INF )
711  ++nlpvars;
712  }
713  assert( nlpvars <= 2 * nlpcons ); /* factor 2 comes from left- and right-hand-sides */
714  }
715 
716  /* compute number of primal scalar variables */
717  nprimalvars = nlpvars + nvarbounds;
718 
719  /* check variable bounds (they should all be nonnegative) */
720  for (i = 0; i < nprimalvars; i++)
721  {
722  viol = MAX(-solvector[i], 0.0);
723  *sumabsviolbnds += viol;
724 
725  if ( viol > *maxabsviolbnds )
726  *maxabsviolbnds = viol;
727 
728 #ifdef SCIP_MORE_DEBUG
729  SCIPdebugMessage("Primal variable %d: val = %g, viol = %g\n", i, solvector[i], viol);
730 #endif
731 
732  if ( solvector[i] < - feastol )
733  {
734  SCIPdebugMessage("primal solution found infeasible (feastol=%f) for primal variable bounds: x[%d] = %g >= 0\n",
735  feastol, i, solvector[i]);
736  *infeasible = TRUE;
737  }
738  }
739 
740  /* compute sizes of primal matrix variables and save maximal size */
741  nprimalmatrixvars = nsdpblocks - nremovedblocks;
742  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &mosekblocksizes, nprimalmatrixvars) );
743  maxblocksize = 0;
744 
745  for (b = 0; b < nsdpblocks; b++)
746  {
747  if ( blockindchanges[b] > -1 )
748  {
749  assert( 0 <= blockindchanges[b] && blockindchanges[b] <= b && (b - blockindchanges[b]) <= (nsdpblocks - nremovedblocks) );
750  mosekblocksizes[b - blockindchanges[b]] = sdpblocksizes[b] - nremovedinds[b];/*lint !e679*/
751 
752  if ( sdpblocksizes[b] - nremovedinds[b] > maxblocksize )
753  maxblocksize = sdpblocksizes[b] - nremovedinds[b];
754  }
755  }
756 
757 
758  /* compute number of primal linear constraints */
759  nprimallpcons = nactivevars;
760 
761  /* check linear constraints */
762  if ( nprimallpcons > 0 )
763  {
764  SCIP_Real* lpconsvals;
765 
766  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &lpconsvals, nprimallpcons) );
767 
768  /* initialize all rows with zero */
769  for (i = 0; i < nprimallpcons; i++)
770  lpconsvals[i] = 0.0;
771 
772  /* first add the matrices A_i */
773  for (b = 0; b < nsdpblocks; b++)
774  {
775  if ( blockindchanges[b] > -1 )
776  {
777  for (blockvar = 0; blockvar < sdpnblockvars[b]; blockvar++)
778  {
779  v = inputtomosekmapper[sdpvar[b][blockvar]];
780 
781  /* check if the variable is active */
782  if ( v > -1 )
783  {
784  assert( v < nactivevars );
785 
786  /* if there are removed indices, we have to adjust the column and row indices accordingly */
787  if ( nremovedinds[b] > 0 )
788  {
789  assert( sdpnblockvarnonz[b][blockvar] <= sdpnnonz );
790  for (k = 0; k < sdpnblockvarnonz[b][blockvar]; k++)
791  {
792  /* rows and cols with active nonzeros should not be removed */
793  assert( 0 <= indchanges[b][sdprow[b][blockvar][k]] && indchanges[b][sdprow[b][blockvar][k]] <= sdprow[b][blockvar][k] );
794  assert( 0 <= indchanges[b][sdpcol[b][blockvar][k]] && indchanges[b][sdpcol[b][blockvar][k]] <= sdpcol[b][blockvar][k] );
795 
796  assert( 0 <= sdprow[b][blockvar][k] && sdprow[b][blockvar][k] < sdpblocksizes[b] );
797  assert( 0 <= sdpcol[b][blockvar][k] && sdpcol[b][blockvar][k] < sdpblocksizes[b] );
798 
799  mosekrow = sdprow[b][blockvar][k] - indchanges[b][sdprow[b][blockvar][k]];
800  mosekcol = sdpcol[b][blockvar][k] - indchanges[b][sdpcol[b][blockvar][k]];
801 
802  /* MOSEK stores matrices in column-first format?!? */
803  blocksize = mosekblocksizes[b - blockindchanges[b]];
804  mosekind = mosekcol * blocksize + mosekrow - mosekcol * (mosekcol + 1) / 2;
805  /* mosekind = mosekrow * (mosekrow + 1) / 2 + mosekcol; */
806 
807  if ( mosekrow == mosekcol )
808  lpconsvals[v] += sdpval[b][blockvar][k] * solmatrices[b - blockindchanges[b]][mosekind];
809  else
810  lpconsvals[v] += 2 * sdpval[b][blockvar][k] * solmatrices[b - blockindchanges[b]][mosekind];
811  }
812  assert( k == sdpnblockvarnonz[b][blockvar] );
813  }
814  else
815  {
816  assert( sdpnblockvarnonz[b][blockvar] <= sdpnnonz );
817  for (k = 0; k < sdpnblockvarnonz[b][blockvar]; k++)
818  {
819  assert( 0 <= sdprow[b][blockvar][k] && sdprow[b][blockvar][k] < sdpblocksizes[b] );
820  assert( 0 <= sdpcol[b][blockvar][k] && sdpcol[b][blockvar][k] < sdpblocksizes[b] );
821 
822  mosekrow = sdprow[b][blockvar][k];
823  mosekcol = sdpcol[b][blockvar][k];
824 
825  /* MOSEK stores matrices in column-first format?!? */
826  blocksize = mosekblocksizes[b - blockindchanges[b]];
827  mosekind = mosekcol * blocksize + mosekrow - mosekcol * (mosekcol + 1) / 2;
828  /* mosekind = mosekrow * (mosekrow + 1) / 2 + mosekcol; */
829 
830  if ( mosekrow == mosekcol )
831  lpconsvals[v] += sdpval[b][blockvar][k] * solmatrices[b - blockindchanges[b]][mosekind];
832  else
833  lpconsvals[v] += 2 * sdpval[b][blockvar][k] * solmatrices[b - blockindchanges[b]][mosekind];
834  }
835  assert( k == sdpnblockvarnonz[b][blockvar] );
836  }
837  }
838  }
839  }
840  }
841 
842  /* add the entries corresponding to the lp-constraints in the dual problem */
843  if ( lpnnonz > 0 )
844  {
845  int currentrow;
846  int varcnt = 0;
847 
848  currentrow = lprow[0];
849  for (nnonz = 0; nnonz < lpnnonz; ++nnonz)
850  {
851  assert( nnonz == 0 || lprow[nnonz-1] <= lprow[nnonz] ); /* rows should be sorted */
852  assert( lprow[nnonz] == currentrow );
853 
854  v = inputtomosekmapper[lpcol[nnonz]];
855  if ( v >= 0 )
856  {
857  assert( v < nactivevars );
858 
859  /* treat left hand side */
860  if ( lplhs[lprow[nnonz]] > - INF )
861  lpconsvals[v] += lpval[nnonz] * solvector[varcnt];
862 
863  /* treat right hand side */
864  if ( lprhs[lprow[nnonz]] < INF )
865  lpconsvals[v] -= lpval[nnonz] * solvector[varcnt + 1];
866  }
867 
868  /* we finished a row */
869  if ( nnonz == lpnnonz - 1 || lprow[nnonz + 1] > currentrow )
870  {
871  /* treat left hand side */
872  if ( lplhs[currentrow] > - INF )
873  varcnt++;
874 
875  /* treat right hand side */
876  if ( lprhs[currentrow] < INF )
877  varcnt++;
878 
879  /* reset counters */
880  if ( nnonz < lpnnonz )
881  currentrow = lprow[nnonz+1];
882  }
883  }
884  assert( varcnt == nlpvars );
885  }
886 
887  /* finally add the entries corresponding to varbounds in the dual problem */
888  for (i = 0; i < nvarbounds; i++)
889  {
890  if ( varboundpos[i] < 0 )
891  {
892  /* lower bound */
893  mosekrow = - varboundpos[i] - 1; /* minus one because we added one to get strictly positive/negative values */
894  assert( 0 <= mosekrow && mosekrow < nactivevars );
895  lpconsvals[mosekrow] += solvector[nlpvars + i];
896  }
897  else
898  {
899  /* upper bound */
900  assert( varboundpos[i] > 0 ); /* we should not have a zero as we wanted a clear differentiation between positive and negative */
901  mosekrow = varboundpos[i] - 1; /* minus one because we added one to get strictly positive/negative values */
902  assert( 0 <= mosekrow && mosekrow < nactivevars );
903  lpconsvals[mosekrow] -= solvector[nlpvars + i];
904  }
905  }
906 
907  /* check all linear constraints for feasibility (they should equal the objective coefficient of the corresponding
908  * dual variable) */
909  for (c = 0; c < nprimallpcons; c++)
910  {
911  viol = REALABS(lpconsvals[c] - objcoefs[c]);
912  *sumabsviolcons += viol;
913 
914  if ( viol > *maxabsviolcons )
915  *maxabsviolcons = viol;
916 
917 #ifdef SCIP_MORE_DEBUG
918  SCIPdebugMessage("Primal linear constraint %d: lhs = %g, rhs = %g, viol = %g\n", c, lpconsvals[c], objcoefs[c], viol);
919 #endif
920 
921  if ( lpconsvals[c] < objcoefs[c] - feastol || lpconsvals[c] > objcoefs[c] + feastol )
922  {
923  SCIPdebugMessage("primal solution found infeasible (feastol=%g) for primal lp constraint: LP-%d = %g == [%g], viol = %g\n",
924  feastol, c, lpconsvals[c], objcoefs[c], viol);
925  *infeasible = TRUE;
926  }
927 
928  }
929  BMSfreeBufferMemoryArray(bufmem, &lpconsvals);
930  }
931 
932  /* check sdp constraints */
933  if ( nprimalmatrixvars > 0 )
934  {
935  SCIP_Real* fullsdpmatrix;
936  SCIP_Real eigenvalue;
937 
938  assert( maxblocksize > 0 );
939 
940  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &fullsdpmatrix, maxblocksize * maxblocksize) );/*lint !e647*/
941 
942  for (b = 0; b < nsdpblocks; b++)
943  {
944  if ( blockindchanges[b] > -1 )
945  {
946  blocksize = mosekblocksizes[b - blockindchanges[b]];
947 
948  /* iterate over all entries in the solution and create corresponding full matrix for LAPACK (in column-first format) */
949  /* NOTE: MOSEK stores matrices in column-first format?!? */
950  for (i = 0; i < blocksize; i++)
951  {
952  for (j = 0; j <= i; j++)
953  {
954  mosekind = j * blocksize + i - j * (j + 1) / 2;
955  fullsdpmatrix[i * blocksize + j] = solmatrices[b - blockindchanges[b]][mosekind];
956  fullsdpmatrix[j * blocksize + i] = solmatrices[b - blockindchanges[b]][mosekind];
957  }
958  }
959 
960  /* compute smallest eigenvalue using LAPACK */
961  SCIP_CALL( SCIPlapackComputeIthEigenvalue(bufmem, FALSE, blocksize, fullsdpmatrix, 1, &eigenvalue, NULL) );
962 
963  viol = MAX(-eigenvalue, 0.0);
964  *sumabsviolsdp += viol;
965 
966  if ( viol > *maxabsviolsdp )
967  *maxabsviolsdp = viol;
968 
969 #ifdef SCIP_MORE_DEBUG
970  SCIPdebugMessage("Primal SDP variable %d: lambda_min = %g, viol = %g\n", b, eigenvalue, viol);
971 #endif
972 
973  if ( eigenvalue < - feastol )
974  {
975  SCIPdebugMessage("primal solution found infeasible (feastol=%.10g) for primal SDP variable %d, smallest eigenvector %.10g\n",
976  feastol, b - blockindchanges[b], eigenvalue);
977  *infeasible = TRUE;
978  }
979  }
980  }
981 
982  BMSfreeBufferMemoryArray(bufmem, &fullsdpmatrix);
983  }
984 
985  BMSfreeBufferMemoryArray(bufmem, &mosekblocksizes);
986  BMSfreeBufferMemoryArray(bufmem, &varboundpos);
987  BMSfreeBufferMemoryArray(bufmem, &objcoefs);
988 
989  return SCIP_OKAY;
990 }
SCIPsdpSolcheckerCheckAndGetViolDual
SCIP_RETCODE SCIPsdpSolcheckerCheckAndGetViolDual(BMS_BUFMEM *bufmem, int nvars, SCIP_Real *lb, SCIP_Real *ub, int nsdpblocks, int *sdpblocksizes, int *sdpnblockvars, int sdpconstnnonz, int *sdpconstnblocknonz, int **sdpconstrow, int **sdpconstcol, SCIP_Real **sdpconstval, int sdpnnonz, int **sdpnblockvarnonz, int **sdpvar, int ***sdprow, int ***sdpcol, SCIP_Real ***sdpval, int **indchanges, int *nremovedinds, int *blockindchanges, int nlpcons, SCIP_Real *lplhs, SCIP_Real *lprhs, int lpnnonz, int *lprow, int *lpcol, SCIP_Real *lpval, SCIP_Real *solvector, SCIP_Real feastol, SCIP_Real epsilon, SCIP_Real *maxabsviolbnds, SCIP_Real *sumabsviolbnds, SCIP_Real *maxabsviolcons, SCIP_Real *sumabsviolcons, SCIP_Real *maxabsviolsdp, SCIP_Real *sumabsviolsdp, SCIP_Bool *infeasible)
Definition: sdpsolchecker.c:281
SCIPsdpSolcheckerCheck
SCIP_RETCODE SCIPsdpSolcheckerCheck(BMS_BUFMEM *bufmem, int nvars, SCIP_Real *lb, SCIP_Real *ub, int nsdpblocks, int *sdpblocksizes, int *sdpnblockvars, int sdpconstnnonz, int *sdpconstnblocknonz, int **sdpconstrow, int **sdpconstcol, SCIP_Real **sdpconstval, int sdpnnonz, int **sdpnblockvarnonz, int **sdpvar, int ***sdprow, int ***sdpcol, SCIP_Real ***sdpval, int **indchanges, int *nremovedinds, int *blockindchanges, int nlpcons, SCIP_Real *lplhs, SCIP_Real *lprhs, int lpnnonz, int *lprow, int *lpcol, SCIP_Real *lpval, SCIP_Real *solvector, SCIP_Real feastol, SCIP_Real epsilon, SCIP_Bool *infeasible)
Definition: sdpsolchecker.c:62
BMS_CALL
#define BMS_CALL(x)
Definition: sdpsolchecker.c:47
sdpsolchecker.h
checks a given SDP solution for feasibility
lapack_interface.h
interface methods for eigenvector computation and matrix multiplication using openblas ...
SCIPlapackComputeIthEigenvalue
SCIP_RETCODE SCIPlapackComputeIthEigenvalue(BMS_BUFMEM *bufmem, SCIP_Bool geteigenvectors, int n, SCIP_Real *A, int i, SCIP_Real *eigenvalue, SCIP_Real *eigenvector)
Definition: lapack_interface.c:174
SCIPsdpSolcheckerCheckAndGetViolPrimal
SCIP_RETCODE SCIPsdpSolcheckerCheckAndGetViolPrimal(BMS_BUFMEM *bufmem, int nvars, SCIP_Real *obj, SCIP_Real *lb, SCIP_Real *ub, int *inputtomosekmapper, int nsdpblocks, int *sdpblocksizes, int *sdpnblockvars, int sdpconstnnonz, int *sdpconstnblocknonz, int **sdpconstrow, int **sdpconstcol, SCIP_Real **sdpconstval, int sdpnnonz, int **sdpnblockvarnonz, int **sdpvar, int ***sdprow, int ***sdpcol, SCIP_Real ***sdpval, int **indchanges, int *nremovedinds, int *blockindchanges, int nremovedblocks, int nlpcons, SCIP_Real *lplhs, SCIP_Real *lprhs, int lpnnonz, int *lprow, int *lpcol, SCIP_Real *lpval, SCIP_Real *solvector, SCIP_Real **solmatrices, SCIP_Real feastol, SCIP_Real epsilon, SCIP_Real *maxabsviolbnds, SCIP_Real *sumabsviolbnds, SCIP_Real *maxabsviolcons, SCIP_Real *sumabsviolcons, SCIP_Real *maxabsviolsdp, SCIP_Real *sumabsviolsdp, SCIP_Bool *infeasible)
Definition: sdpsolchecker.c:550
INF
#define INF
Definition: sdpsolchecker.c:57
Generated by   doxygen 1.8.11