SCIP-SDP  3.2.0
cons_sdp.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-2020 Discrete Optimization, TU Darmstadt */
9 /* */
10 /* */
11 /* This program is free software; you can redistribute it and/or */
12 /* modify it under the terms of the GNU Lesser General Public License */
13 /* as published by the Free Software Foundation; either version 3 */
14 /* of the License, or (at your option) any later version. */
15 /* */
16 /* This program is distributed in the hope that it will be useful, */
17 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
18 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
19 /* GNU Lesser General Public License for more details. */
20 /* */
21 /* You should have received a copy of the GNU Lesser General Public License */
22 /* along with this program; if not, write to the Free Software */
23 /* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.*/
24 /* */
25 /* */
26 /* Based on SCIP - Solving Constraint Integer Programs */
27 /* Copyright (C) 2002-2020 Zuse Institute Berlin */
28 /* SCIP is distributed under the terms of the SCIP Academic Licence, */
29 /* see file COPYING in the SCIP distribution. */
30 /* */
31 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
32 
49 /* #define SCIP_DEBUG */
50 /* #define SCIP_MORE_DEBUG /\* shows all cuts added and prints constraint after parsing *\/ */
51 /* #define PRINT_HUMAN_READABLE /\* change the output of PRINTCONS to a better readable format (dense instead of sparse), WHICH CAN NO LONGER BE PARSED *\/ */
52 /* #define PRINTMATRICES /\* Should all matrices appearing in best rank-1 approximation heuristic be printed? *\/ */
53 
54 #include "cons_sdp.h"
55 
56 #include <assert.h> /*lint !e451*/
57 #include <string.h> /* for NULL, strcmp */
58 #include <ctype.h> /* for isspace */
59 #include <math.h>
60 #include "sdpi/lapack_interface.h"
61 
62 #include "scipsdp/SdpVarmapper.h"
63 #include "scipsdp/SdpVarfixer.h"
64 
65 #include "scip/cons_linear.h" /* for SCIPcreateConsLinear */
66 #include "scip/cons_quadratic.h" /* for SCIPcreateConsBasicQuadratic */
67 #include "scip/scip_cons.h" /* for SCIPgetConsVars */
68 #include "scip/scip.h" /* for SCIPallocBufferArray, etc */
69 #include "scip/def.h"
70 
71 #ifdef OMP
72 #include "omp.h" /* for changing the number of threads */
73 #endif
74 
75 /* turn off lint warnings for whole file: */
76 /*lint --e{788,818}*/
77 
78 #define CONSHDLR_NAME "SDP"
79 #define CONSHDLR_DESC "SDP constraints of the form \\sum_{j} A_j y_j - A_0 psd"
80 
81 #define CONSHDLRRANK1_NAME "SDPrank1"
82 #define CONSHDLRRANK1_DESC "rank 1 SDP constraints"
83 
84 #define CONSHDLR_SEPAPRIORITY +1000000
85 #define CONSHDLR_ENFOPRIORITY -2000000
86 #define CONSHDLR_CHECKPRIORITY -2000000
87 #define CONSHDLR_SEPAFREQ 1
88 #define CONSHDLR_EAGERFREQ 100
90 #define CONSHDLR_MAXPREROUNDS -1
91 #define CONSHDLR_DELAYSEPA FALSE
92 #define CONSHDLR_NEEDSCONS TRUE
94 #define CONSHDLR_PRESOLTIMING SCIP_PRESOLTIMING_FAST
95 #define PARSE_STARTSIZE 1
96 #define PARSE_SIZEFACTOR 10
97 #define DEFAULT_DIAGGEZEROCUTS TRUE
98 #define DEFAULT_DIAGZEROIMPLCUTS TRUE
99 #define DEFAULT_TWOMINORLINCONSS FALSE
100 #define DEFAULT_TWOMINORPRODCONSS FALSE
101 #define DEFAULT_QUADCONSRANK1 TRUE
102 #define DEFAULT_UPGRADQUADCONSS FALSE
103 #define DEFAULT_UPGRADEKEEPQUAD FALSE
104 #define DEFAULT_MAXNVARSQUADUPGD 1000
105 #define DEFAULT_RANK1APPROXHEUR FALSE
106 #ifdef OMP
107 #define DEFAULT_NTHREADS 1
108 #endif
109 
111 struct SCIP_ConsData
112 {
113  int nvars;
114  int nnonz;
115  int blocksize;
116  int* nvarnonz;
117  int** col;
118  int** row;
119  SCIP_Real** val;
120  SCIP_VAR** vars;
121  int* locks;
122  int constnnonz;
123  int* constcol;
124  int* constrow;
125  SCIP_Real* constval;
126  SCIP_Real maxrhsentry;
127  SCIP_Bool rankone;
128  int* maxevsubmat;
129  SCIP_Bool addedquadcons;
130 };
131 
133 struct SCIP_ConshdlrData
134 {
135  int neigveccuts;
136  SCIP_Bool diaggezerocuts;
137  int ndiaggezerocuts;
138  int n1x1blocks;
139  SCIP_Bool diagzeroimplcuts;
140  SCIP_Bool twominorlinconss;
141  SCIP_Bool twominorprodconss;
142  SCIP_Bool quadconsrank1;
143  SCIP_Bool upgradquadconss;
144  SCIP_Bool upgradekeepquad;
145  int maxnvarsquadupgd;
146  SCIP_Bool triedlinearconss;
147  SCIP_Bool rank1approxheur;
148 #ifdef OMP
149  int nthreads;
150 #endif
151  int* quadconsidx;
152  SCIP_VAR** quadconsvars;
153  int nquadconsidx;
154  SCIP_VAR*** X;
155  int nsdpvars;
156  SCIP_CONS* sdpcons;
157  SCIP_CONSHDLRDATA* sdpconshdlrdata;
158 };
159 
163 static
164 SCIP_RETCODE convertRowToColFormatFullMatrix(
165  int rows,
166  int cols,
167  SCIP_Real* rowmatrix,
168  SCIP_Real* colmatrix
170  )
171 {
172  int i;
173  int j;
174  SCIP_Real act;
175 
176  assert( rows > 0 );
177  assert( cols > 0 );
178  assert( rowmatrix != NULL );
179  assert( colmatrix != NULL );
180 
181  for (i = 0; i < rows; ++i)
182  {
183  for (j = 0; j < cols; ++j)
184  {
185  act = rowmatrix[i*cols + j];
186  colmatrix[j*rows + i] = act;
187  }
188  }
189 
190  return SCIP_OKAY;
191 }
192 
194 static
195 SCIP_RETCODE scaleRowsMatrix(
196  int blocksize, /* number of rows and columns */
197  SCIP_Real* matrix, /* matrix entries given as blocksize^2 array */
198  SCIP_Real* scale /* array of length blocksize to multiply the rows of matrix with */
199  )
200 {
201  int r; /* rows */
202  int c; /* columns */
203 
204  assert( blocksize >= 0 );
205  assert( matrix != NULL );
206  assert( scale != NULL );
207 
208  for (r = 0; r < blocksize; r++)
209  {
210  for (c = 0; c < blocksize; c++)
211  {
212  /* row-first format! */
213  matrix[r * blocksize + c] *= scale[r]; /*lint !e679*/
214  }
215  }
216 
217  return SCIP_OKAY;
218 }
219 
221 static
222 SCIP_RETCODE expandSymMatrix(
223  int size,
224  SCIP_Real* symMat,
225  SCIP_Real* fullMat
226  )
227 {
228  int i;
229  int j;
230  int ind = 0;
231 
232  assert( size >= 0 );
233  assert( symMat != NULL );
234  assert( fullMat != NULL );
235 
236  /* traverse the lower triangular part in the order of the indices and copy the values to both lower and upper triangular part */
237  for (i = 0; i < size; i++)
238  {
239  for (j = 0; j <= i; j++)
240  {
241  assert( ind == SCIPconsSdpCompLowerTriangPos(i,j) );
242  fullMat[i*size + j] = symMat[ind]; /*lint !e679*/
243  fullMat[j*size + i] = symMat[ind]; /*lint !e679*/
244  ind++;
245  }
246  }
247 
248  return SCIP_OKAY;
249 }
250 
254 static
255 SCIP_RETCODE computeSdpMatrix(
256  SCIP* scip,
257  SCIP_CONS* cons,
258  SCIP_SOL* y,
259  SCIP_Real* matrix
260  )
261 {
262  SCIP_CONSDATA* consdata;
263  SCIP_Real yval;
264  int blocksize;
265  int nvars;
266  int ind;
267  int i;
268 
269  assert( cons != NULL );
270  assert( matrix != NULL );
271 
272  consdata = SCIPconsGetData(cons);
273  nvars = consdata->nvars;
274  blocksize = consdata->blocksize;
275 
276  /* initialize the matrix with 0 */
277  for (i = 0; i < (blocksize * (blocksize + 1))/2; i++)
278  matrix[i] = 0.0;
279 
280  /* add the non-constant-part */
281  for (i = 0; i < nvars; i++)
282  {
283  yval = SCIPgetSolVal(scip, y, consdata->vars[i]);
284  for (ind = 0; ind < consdata->nvarnonz[i]; ind++)
285  matrix[SCIPconsSdpCompLowerTriangPos(consdata->row[i][ind], consdata->col[i][ind])] += yval * consdata->val[i][ind];
286  }
287 
288  /* substract the constant part */
289  for (ind = 0; ind < consdata->constnnonz; ind++)
290  matrix[SCIPconsSdpCompLowerTriangPos(consdata->constrow[ind], consdata->constcol[ind])] -= consdata->constval[ind];
291 
292  return SCIP_OKAY;
293 }
294 
296 static
297 SCIP_RETCODE isMatrixRankOne(
298  SCIP* scip,
299  SCIP_CONS* cons,
300  SCIP_SOL* sol,
301  SCIP_Bool* result
302  )
303 {
304  SCIP_CONSDATA* consdata;
305  SCIP_Real* matrix = NULL;
306  SCIP_Real* fullmatrix = NULL;
307  SCIP_Real eigenvalue;
308  int blocksize;
309  SCIP_RESULT resultSDPtest;
310  int i;
311  int j;
312  int ind1 = 0;
313  int ind2 = 0;
314  SCIP_Real submatrix[4];
315  SCIP_Real largestminev = 0.0;
316 
317  assert( cons != NULL );
318  assert( result != NULL );
319 
320  consdata = SCIPconsGetData(cons);
321  assert( consdata != NULL );
322 
323  blocksize = consdata->blocksize;
324 
325  resultSDPtest = SCIP_INFEASIBLE;
326 
327  SCIP_CALL( SCIPconsSdpCheckSdpCons(scip, cons, sol, FALSE, &resultSDPtest) );
328 
329  if ( resultSDPtest == SCIP_INFEASIBLE )
330  {
331  SCIPerrorMessage("Try to check for a matrix to be rank 1 even if the matrix is not psd.\n");
332  return SCIP_ERROR;
333  }
334 
335  /* allocate memory to store full matrix */
336  SCIP_CALL( SCIPallocBufferArray(scip, &matrix, (blocksize * (blocksize+1))/2 ) );
337  SCIP_CALL( SCIPallocBufferArray(scip, &fullmatrix, blocksize * blocksize ) );
338 
339  /* compute the matrix \f$ \sum_j A_j y_j - A_0 \f$ */
340  SCIP_CALL( computeSdpMatrix(scip, cons, sol, matrix) );
341 
342  /* expand it because LAPACK wants the full matrix instead of the lower triangular part */
343  SCIP_CALL( expandSymMatrix(blocksize, matrix, fullmatrix) );
344 
345  /* compute the second largest eigenvalue */
346  SCIP_CALL( SCIPlapackComputeIthEigenvalue(SCIPbuffer(scip), FALSE, blocksize, fullmatrix, blocksize - 1, &eigenvalue, NULL) );
347 
348  /* the matrix is rank 1 iff the second largest eigenvalue is zero (since the matrix is symmetric and psd) */
349 
350  /* changed eigenvalue to 0.1*eigenvalue here, since here seems to be a problem with optimal solutions not satisfying
351  SCIPisFeasEQ(scip, 0.1*eigenvalue, 0.0), even if they are feasible for all other constraints, including the
352  quadratic 2x2 principal minor constraints. */
353  if ( SCIPisFeasEQ(scip, 0.1*eigenvalue, 0.0) )
354  *result = TRUE;
355  else
356  {
357  *result = FALSE;
358 
359  /* if the matrix is not rank 1, compute minimal eigenvalues of 2x2 minors */
360  for (i = 0; i < blocksize; ++i)
361  {
362  for (j = 0; j < i; ++j)
363  {
364  submatrix[0] = matrix[SCIPconsSdpCompLowerTriangPos(i,i)];
365  submatrix[1] = matrix[SCIPconsSdpCompLowerTriangPos(i,j)];
366  submatrix[2] = matrix[SCIPconsSdpCompLowerTriangPos(i,j)];
367  submatrix[3] = matrix[SCIPconsSdpCompLowerTriangPos(j,j)];
368 
369  SCIP_CALL( SCIPlapackComputeIthEigenvalue(SCIPbuffer(scip), FALSE, 2, submatrix, 1, &eigenvalue, NULL) );
370  /* TODO: Compute eigenvalues by solving quadratic constraint */
371 
372  if ( eigenvalue > largestminev )
373  {
374  largestminev = eigenvalue;
375  ind1 = i;
376  ind2 = j;
377  }
378  }
379  }
380 
381  assert( ind1 > 0 || ind2 > 0 );
382  assert( SCIPisFeasPositive(scip, largestminev) ); /* because the second largest eigenvalue satisfies SCIPisFeasPositive. */
383 
384  /* save indices for submatrix with largest minimal eigenvalue */
385  consdata->maxevsubmat[0] = ind1;
386  consdata->maxevsubmat[1] = ind2;
387  }
388 
389  SCIPfreeBufferArray(scip, &fullmatrix);
390  SCIPfreeBufferArray(scip, &matrix);
391 
392  return SCIP_OKAY;
393 }
394 
396 static
397 SCIP_RETCODE multiplyConstraintMatrix(
398  SCIP_CONS* cons,
399  int j,
400  SCIP_Real* v,
401  SCIP_Real* vAv
402  )
403 {
404  SCIP_CONSDATA* consdata;
405  int i;
406 
407  assert( cons != NULL );
408  assert( j >= 0 );
409  assert( vAv != NULL );
410 
411  consdata = SCIPconsGetData(cons);
412  assert( consdata != NULL );
413 
414  assert( j < consdata->nvars );
415 
416  /* initialize the product with 0 */
417  *vAv = 0.0;
418 
419  for (i = 0; i < consdata->nvarnonz[j]; i++)
420  {
421  if (consdata->col[j][i] == consdata->row[j][i])
422  *vAv += v[consdata->col[j][i]] * consdata->val[j][i] * v[consdata->row[j][i]];
423  else
424  {
425  /* Multiply by 2, because the matrix is symmetric and there is one identical contribution each from lower and upper triangular part. */
426  *vAv += 2.0 * v[consdata->col[j][i]] * consdata->val[j][i] * v[consdata->row[j][i]];
427  }
428  }
429 
430  return SCIP_OKAY;
431 }
432 
436 static
437 SCIP_RETCODE setMaxRhsEntry(
438  SCIP_CONS* cons
439  )
440 {
441  SCIP_CONSDATA* consdata;
442  SCIP_Real max;
443  int i;
444 
445  consdata = SCIPconsGetData(cons);
446  assert( consdata != NULL );
447 
448  /* initialize max with zero (this is used if there is no constant-matrix) */
449  max = 0.0;
450 
451  /* iterate over the entries of the constant matrix, updating max if a higher absolute value is found */
452  for (i = 0; i < consdata->constnnonz; i++)
453  {
454  if ( REALABS(consdata->constval[i]) > max )
455  max = REALABS(consdata->constval[i]);
456  }
457 
458  consdata->maxrhsentry = max;
459 
460  return SCIP_OKAY;
461 }
462 
463 
469 static
470 SCIP_RETCODE cutUsingEigenvector(
471  SCIP* scip,
472  SCIP_CONS* cons,
473  SCIP_SOL* sol,
474  SCIP_Real* coeff,
475  SCIP_Real* lhs
476  )
477 {
478  SCIP_CONSDATA* consdata;
479  SCIP_Real* eigenvalues = NULL;
480  SCIP_Real* matrix = NULL;
481  SCIP_Real* fullmatrix = NULL;
482  SCIP_Real* fullconstmatrix = NULL;
483  SCIP_Real* eigenvector = NULL;
484  SCIP_Real* output_vector = NULL;
485  int blocksize;
486  int j;
487 
488  assert( cons != NULL );
489  assert( lhs != NULL );
490 
491  *lhs = 0.0;
492 
493  consdata = SCIPconsGetData(cons);
494  assert( consdata != NULL );
495 
496  blocksize = consdata->blocksize;
497 
498  SCIP_CALL( SCIPallocBufferArray(scip, &eigenvalues, blocksize) );
499  SCIP_CALL( SCIPallocBufferArray(scip, &matrix, (blocksize * (blocksize+1))/2 ) ); /*lint !e647*/
500  SCIP_CALL( SCIPallocBufferArray(scip, &fullmatrix, blocksize * blocksize ) ); /*lint !e647*/
501  SCIP_CALL( SCIPallocBufferArray(scip, &fullconstmatrix, blocksize * blocksize) ); /*lint !e647*/
502  SCIP_CALL( SCIPallocBufferArray(scip, &eigenvector, blocksize) );
503  SCIP_CALL( SCIPallocBufferArray(scip, &output_vector, blocksize) );
504 
505  /* compute the matrix \f$ \sum_j A_j y_j - A_0 \f$ */
506  SCIP_CALL( computeSdpMatrix(scip, cons, sol, matrix) );
507 
508  /* expand it because LAPACK wants the full matrix instead of the lower triangular part */
509  SCIP_CALL( expandSymMatrix(blocksize, matrix, fullmatrix) );
510 
511  SCIP_CALL( SCIPlapackComputeIthEigenvalue(SCIPbuffer(scip), TRUE, blocksize, fullmatrix, 1, eigenvalues, eigenvector) );
512 
513  /* get full constant matrix */
514  SCIP_CALL( SCIPconsSdpGetFullConstMatrix(scip, cons, fullconstmatrix) );
515 
516  /* multiply eigenvector with constant matrix to get lhs (after multiplying again with eigenvector from the left) */
517  SCIP_CALL( SCIPlapackMatrixVectorMult(blocksize, blocksize, fullconstmatrix, eigenvector, output_vector) );
518 
519  for (j = 0; j < blocksize; ++j)
520  *lhs += eigenvector[j] * output_vector[j];
521 
522  /* compute \f$ v^T A_j v \f$ for eigenvector v and each matrix \f$ A_j \f$ to get the coefficients of the LP cut */
523  for (j = 0; j < consdata->nvars; ++j)
524  {
525  SCIP_CALL( multiplyConstraintMatrix(cons, j, eigenvector, &coeff[j]) );
526  }
527 
528  SCIPfreeBufferArray(scip, &output_vector);
529  SCIPfreeBufferArray(scip, &eigenvector);
530  SCIPfreeBufferArray(scip, &fullconstmatrix);
531  SCIPfreeBufferArray(scip, &fullmatrix);
532  SCIPfreeBufferArray(scip, &matrix);
533  SCIPfreeBufferArray(scip, &eigenvalues);
534 
535  return SCIP_OKAY;
536 }
537 
539 SCIP_RETCODE SCIPconsSdpCheckSdpCons(
540  SCIP* scip,
541  SCIP_CONS* cons,
542  SCIP_SOL* sol,
543  SCIP_Bool printreason,
544  SCIP_RESULT* result
545  )
546 { /*lint --e{715}*/
547  SCIP_CONSDATA* consdata;
548  int blocksize;
549  SCIP_Real check_value;
550  SCIP_Real eigenvalue;
551  SCIP_Real* matrix = NULL;
552  SCIP_Real* fullmatrix = NULL;
553 
554  assert( scip != NULL );
555  assert( cons != NULL );
556  assert( result != NULL );
557 
558  consdata = SCIPconsGetData(cons);
559  assert( consdata != NULL );
560  assert( consdata->rankone || strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(cons)), CONSHDLR_NAME) == 0 );
561  assert( ! consdata->rankone || strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(cons)), CONSHDLRRANK1_NAME) == 0 );
562  blocksize = consdata->blocksize;
563 
564  SCIP_CALL( SCIPallocBufferArray(scip, &matrix, (blocksize * (blocksize+1)) / 2) ); /*lint !e647*/
565  SCIP_CALL( SCIPallocBufferArray(scip, &fullmatrix, blocksize * blocksize) ); /*lint !e647*/
566  SCIP_CALL( computeSdpMatrix(scip, cons, sol, matrix) );
567  SCIP_CALL( expandSymMatrix(blocksize, matrix, fullmatrix) );
568 
569  SCIP_CALL( SCIPlapackComputeIthEigenvalue(SCIPbuffer(scip), FALSE, blocksize, fullmatrix, 1, &eigenvalue, NULL) );
570 
571  /* This enables checking the second DIMACS Error Norm: err=max{0, -lambda_min(x)/(1+maximumentry of rhs)}, in that case it also needs
572  * to be changed in the sdpi (and implemented there first), when checking feasibility of problems where all variables are fixed */
573 #ifdef DIMACS
574  check_value = (-eigenvalue) / (1.0 + consdata->maxrhsentry);
575 #else
576  check_value = -eigenvalue;
577 #endif
578 
579  if ( SCIPisFeasLE(scip, check_value, 0.0) )
580  *result = SCIP_FEASIBLE;
581  else
582  {
583  *result = SCIP_INFEASIBLE;
584  if ( printreason )
585  {
586  SCIPinfoMessage(scip, NULL, "SDP-constraint <%s> violated: non psd matrix (eigenvalue %f, dimacs error norm = %f).\n", SCIPconsGetName(cons), eigenvalue, check_value);
587  SCIP_CALL( SCIPprintCons(scip, cons, NULL) );
588  }
589  }
590 
591  if ( sol != NULL )
592  SCIPupdateSolConsViolation(scip, sol, -eigenvalue, (-eigenvalue) / (1.0 + consdata->maxrhsentry));
593 
594  SCIPfreeBufferArray(scip, &fullmatrix);
595  SCIPfreeBufferArray(scip, &matrix);
596 
597  return SCIP_OKAY;
598 }
599 
601 static
602 SCIP_RETCODE separateSol(
603  SCIP* scip,
604  SCIP_CONSHDLR* conshdlr,
605  SCIP_CONS* cons,
606  SCIP_SOL* sol,
607  SCIP_RESULT* result
608  )
609 {
610  char cutname[SCIP_MAXSTRLEN];
611  SCIP_CONSDATA* consdata;
612  SCIP_CONSHDLRDATA* conshdlrdata;
613  SCIP_Real lhs = 0.0;
614  SCIP_Real* coeff = NULL;
615  SCIP_COL** cols;
616  SCIP_Real* vals;
617  SCIP_ROW* row;
618  int nvars;
619  int j;
620  int len;
621 #ifndef NDEBUG
622  int snprintfreturn; /* this is used to assert that the SCIP string concatenation works */
623 #endif
624 
625  assert( cons != NULL );
626 
627  consdata = SCIPconsGetData(cons);
628  assert( consdata != NULL );
629 
630  nvars = consdata->nvars;
631  SCIP_CALL( SCIPallocBufferArray(scip, &coeff, nvars ) );
632 
633  SCIP_CALL( cutUsingEigenvector(scip, cons, sol, coeff, &lhs) );
634 
635  SCIP_CALL( SCIPallocBufferArray(scip, &cols, nvars ) );
636  SCIP_CALL( SCIPallocBufferArray(scip, &vals, nvars ) );
637 
638  len = 0;
639  for (j = 0; j < nvars; ++j)
640  {
641  if ( SCIPisZero(scip, coeff[j]) )
642  continue;
643 
644  cols[len] = SCIPvarGetCol(SCIPvarGetProbvar(consdata->vars[j]));
645  vals[len] = coeff[j];
646  ++len;
647  }
648 
649  conshdlrdata = SCIPconshdlrGetData(conshdlr);
650  assert( conshdlrdata != NULL );
651 
652 #ifndef NDEBUG
653  snprintfreturn = SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "sepa_eig_sdp_%d", ++(conshdlrdata->neigveccuts));
654  assert( snprintfreturn < SCIP_MAXSTRLEN ); /* check whether name fit into string */
655 #else
656  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "sepa_eig_sdp_%d", ++(conshdlrdata->neigveccuts));
657 #endif
658 
659 #if ( SCIP_VERSION >= 700 || (SCIP_VERSION >= 602 && SCIP_SUBVERSION > 0) )
660  SCIP_CALL( SCIPcreateRowConshdlr(scip, &row, conshdlr, cutname , len, cols, vals, lhs, SCIPinfinity(scip), FALSE, FALSE, TRUE) );
661 #else
662  SCIP_CALL( SCIPcreateRowCons(scip, &row, conshdlr, cutname , len, cols, vals, lhs, SCIPinfinity(scip), FALSE, FALSE, TRUE) );
663 #endif
664 
665  if ( SCIPisCutEfficacious(scip, sol, row) )
666  {
667  SCIP_Bool infeasible;
668 #ifdef SCIP_MORE_DEBUG
669  SCIPinfoMessage(scip, NULL, "Added cut %s: ", cutname);
670  SCIPinfoMessage(scip, NULL, "%f <= ", lhs);
671  for (j = 0; j < len; j++)
672  SCIPinfoMessage(scip, NULL, "+ (%f)*%s", vals[j], SCIPvarGetName(SCIPcolGetVar(cols[j])));
673  SCIPinfoMessage(scip, NULL, "\n");
674 #endif
675 
676  SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
677  if ( infeasible )
678  *result = SCIP_CUTOFF;
679  else
680  *result = SCIP_SEPARATED;
681  SCIP_CALL( SCIPresetConsAge(scip, cons) );
682  }
683 
684  SCIP_CALL( SCIPreleaseRow(scip, &row) );
685 
686  SCIPfreeBufferArray(scip, &vals);
687  SCIPfreeBufferArray(scip, &cols);
688  SCIPfreeBufferArray(scip, &coeff);
689 
690  return SCIP_OKAY;
691 }
692 
696 static
697 SCIP_RETCODE diagGEzero(
698  SCIP* scip,
699  SCIP_CONSHDLR* conshdlr,
700  SCIP_CONS** conss,
701  int nconss,
702  int* naddconss,
703  int* nchgbds,
704  SCIP_RESULT* result
705  )
706 {
707  char cutname[SCIP_MAXSTRLEN];
708  SCIP_CONSHDLRDATA* conshdlrdata;
709  SCIP_CONSDATA* consdata;
710  SCIP_Real* matrix;
711  SCIP_Real* consvals;
712  SCIP_VAR** consvars;
713  SCIP_Real* diagentries;
714  int blocksize;
715  int nvars;
716  int i;
717  int j;
718  int k;
719  int c;
720 
721  assert( scip != NULL );
722  assert( naddconss != NULL );
723  assert( nchgbds != NULL );
724  assert( result != NULL );
725  assert( *result != SCIP_CUTOFF );
726 
727  conshdlrdata = SCIPconshdlrGetData(conshdlr);
728  assert( conshdlrdata != NULL );
729 
730  for (c = 0; c < nconss; ++c)
731  {
732  consdata = SCIPconsGetData(conss[c]);
733  assert( consdata != NULL );
734  assert( consdata->rankone || strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(conss[c])), CONSHDLR_NAME) == 0 );
735  assert( ! consdata->rankone || strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(conss[c])), CONSHDLRRANK1_NAME) == 0 );
736 
737  blocksize = consdata->blocksize;
738  nvars = consdata->nvars;
739 
740  SCIP_CALL( SCIPallocBufferArray(scip, &consvars, nvars) );
741  SCIP_CALL( SCIPallocBufferArray(scip, &consvals, nvars) );
742 
743  SCIP_CALL( SCIPallocBufferArray(scip, &matrix, (blocksize * (blocksize + 1)) / 2) ); /*lint !e647*/
744  SCIP_CALL( SCIPconsSdpGetLowerTriangConstMatrix(scip, conss[c], matrix) );
745 
746  /* allocate diagonal entries and intit to 0.0 */
747  SCIP_CALL( SCIPallocClearBufferArray(scip, &diagentries, blocksize * nvars) ); /*lint !e647*/
748 
749  /* get the (k,k)-entry of every matrix A_j */
750  for (j = 0; j < nvars; ++j)
751  {
752  /* go through the nonzeros of A_j and look for diagonal entries */
753  for (i = 0; i < consdata->nvarnonz[j]; i++)
754  {
755  if ( consdata->col[j][i] == consdata->row[j][i] )
756  diagentries[consdata->col[j][i] * nvars + j] = consdata->val[j][i]; /*lint !e679*/
757  }
758  }
759 
760  /* add the LP-cuts to SCIP */
761  for (k = 0; k < blocksize && *result != SCIP_CUTOFF; ++k)
762  {
763  SCIP_CONS* cons;
764  SCIP_Real lhs;
765  SCIP_Real activitylb = 0.0;
766  int cnt = 0;
767 
768  for ( j = 0; j < nvars; ++j)
769  {
770  SCIP_Real val;
771  SCIP_VAR* var;
772 
773  val = diagentries[k * nvars + j];
774  if ( ! SCIPisZero(scip, val) )
775  {
776  var = consdata->vars[j];
777  consvals[cnt] = val;
778  consvars[cnt++] = var;
779  /* compute lower bound on activity */
780  if ( val > 0 )
781  activitylb += val * SCIPvarGetLbGlobal(var);
782  else
783  activitylb += val * SCIPvarGetUbGlobal(var);
784  }
785  }
786 
787  /* the lhs is the (k,k)-entry of the constant matrix */
788  lhs = matrix[SCIPconsSdpCompLowerTriangPos(k,k)];
789 
790  if ( cnt == 0 )
791  {
792  /* if there are no variables, but the lhs is positive, we are infeasible */
793  if ( SCIPisPositive(scip, lhs) )
794  *result = SCIP_CUTOFF;
795  }
796  else if ( cnt == 1 )
797  {
798  SCIP_Bool infeasible = FALSE;
799  SCIP_Bool tightened;
800  SCIP_VAR* var;
801  SCIP_Real val;
802 
803  var = consvars[0];
804  val = consvals[0];
805  assert( var != NULL );
806 
807  /* try to tighten bound */
808  if ( SCIPisPositive(scip, val) )
809  {
810  SCIP_CALL( SCIPtightenVarLb(scip, var, lhs / val, FALSE, &infeasible, &tightened) );
811  if ( tightened )
812  {
813  SCIPdebugMsg(scip, "Tightend lower bound of <%s> to %g because of diagonal values of SDP-constraint %s!\n",
814  SCIPvarGetName(var), SCIPvarGetLbGlobal(var), SCIPconsGetName(conss[c]));
815  ++(*nchgbds);
816  }
817  }
818  else if ( SCIPisNegative(scip, val) )
819  {
820  SCIP_CALL( SCIPtightenVarUb(scip, var, lhs / val, FALSE, &infeasible, &tightened) );
821  if ( tightened )
822  {
823  SCIPdebugMsg(scip, "Tightend upper bound of <%s> to %g because of diagonal values of SDP-constraint %s!\n",
824  SCIPvarGetName(var), SCIPvarGetUbGlobal(var), SCIPconsGetName(conss[c]));
825  ++(*nchgbds);
826  }
827  }
828 
829  if ( infeasible )
830  *result = SCIP_CUTOFF;
831  }
832  /* generate linear inequality if lower bound on activity is less than the lhs, so the cut is not redundant */
833  else if ( SCIPisLT(scip, activitylb, lhs) )
834  {
835  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "diag_ge_zero_%d", ++(conshdlrdata->ndiaggezerocuts));
836 
837  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, cutname, cnt, consvars, consvals, lhs, SCIPinfinity(scip),
838  TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE) ); /*lint !e679*/
839 
840  SCIP_CALL( SCIPaddCons(scip, cons) );
841 #ifdef SCIP_MORE_DEBUG
842  SCIPinfoMessage(scip, NULL, "Added diagonal lp-constraint: ");
843  SCIP_CALL( SCIPprintCons(scip, cons, NULL) );
844  SCIPinfoMessage(scip, NULL, "\n");
845 #endif
846  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
847  ++(*naddconss);
848  }
849  }
850 
851  SCIPfreeBufferArray(scip, &diagentries);
852  SCIPfreeBufferArray(scip, &matrix);
853  SCIPfreeBufferArray(scip, &consvals);
854  SCIPfreeBufferArray(scip, &consvars);
855  }
856 
857  return SCIP_OKAY;
858 }
859 
867 static
868 SCIP_RETCODE diagZeroImpl(
869  SCIP* scip,
870  SCIP_CONS** conss,
871  int nconss,
872  int* naddconss
873  )
874 {
875  char cutname[SCIP_MAXSTRLEN];
876  /* if entry k is >= 0, gives the number of non-diagonal nonzero-entries in row k of the constant matrix (and the number
877  * of entries in constnonzeroentries), if -1 C_kk =!= 0 (which means we didnot allocate memory for diagvars), if -2
878  * either A_jk =!= 0 for all j with C_jk =!= 0 or A_kk =!= 0 for some continuous variable (which means we did allocate
879  * memory for diagvars but cannot use the cut */
880  int* nconstnonzeroentries;
881  int** constnonzeroentries;
882  SCIP_CONSDATA* consdata;
883  SCIP_CONS* cons;
884  SCIP_VAR** vars;
885  SCIP_Real* vals;
886  int** diagvars;
887  int* ndiagvars;
888  int blocksize;
889  int i;
890  int j;
891  int nvars;
892  int v;
893  int k;
894  int l;
895  SCIP_Bool anycutvalid;
896 #ifndef NDEBUG
897  int snprintfreturn;
898 #endif
899 
900  assert( scip != NULL );
901  assert( naddconss != NULL );
902 
903  for (i = 0; i < nconss; ++i)
904  {
905  assert( conss[i] != NULL );
906  consdata = SCIPconsGetData(conss[i]);
907  assert( consdata != NULL );
908 
909  blocksize = consdata->blocksize;
910  nvars = consdata->nvars;
911  SCIP_CALL( SCIPallocBufferArray(scip, &nconstnonzeroentries, blocksize) );
912  SCIP_CALL( SCIPallocBufferArray(scip, &constnonzeroentries, blocksize) );
913  for (j = 0; j < blocksize; j++)
914  {
915  nconstnonzeroentries[j] = 0;
916  SCIP_CALL( SCIPallocBufferArray(scip, &constnonzeroentries[j], blocksize) );
917  }
918 
919  /* iterate over all nonzeros of the constant matrix and check which diagonal and non-diagonal entries are nonzero */
920  for (j = 0; j < consdata->constnnonz; j++)
921  {
922  /* if it is a nondiagonal-entry we add this row/column to the constnonzeroentries entries unless we already found a
923  * diagonal entry for this row/column */
924  if ( (consdata->constcol[j] != consdata->constrow[j]) )
925  {
926  assert( ! SCIPisZero(scip, consdata->constval[j]) );
927  if ( nconstnonzeroentries[consdata->constcol[j]] >= 0 )
928  {
929  constnonzeroentries[consdata->constcol[j]][nconstnonzeroentries[consdata->constcol[j]]] = consdata->constrow[j];
930  nconstnonzeroentries[consdata->constcol[j]]++;
931  }
932 
933  if ( nconstnonzeroentries[consdata->constrow[j]] >= 0 )
934  {
935  constnonzeroentries[consdata->constrow[j]][nconstnonzeroentries[consdata->constrow[j]]] = consdata->constcol[j];
936  nconstnonzeroentries[consdata->constrow[j]]++;
937  }
938  }
939  /* if we find a diagonal entry in the constant matrix, we remember that we cannot add a cut for this index */
940  else
941  {
942  assert( ! SCIPisZero(scip, consdata->constval[j]) );
943  nconstnonzeroentries[consdata->constcol[j]] = -1;
944  }
945  }
946 
947  /* diagvars[j] is an array with all variables with a diagonal entry (j,j) in the corresponding matrix, if nconstnonzeroentries[j] =!= -1 or NULL otherwise
948  * the outer array goes over all rows to ease the access, but only for those that are really needed memory will be allocated */
949  SCIP_CALL( SCIPallocBufferArray(scip, &diagvars, blocksize) );
950  SCIP_CALL( SCIPallocBufferArray(scip, &ndiagvars, blocksize) );
951  anycutvalid = FALSE;
952  for (j = 0; j < blocksize; ++j)
953  {
954  ndiagvars[j] = 0;
955  if ( nconstnonzeroentries[j] > 0 )
956  {
957  SCIP_CALL( SCIPallocBufferArray(scip, &(diagvars[j]), nvars) );
958  anycutvalid = TRUE;
959  }
960  }
961 
962  /* if no cuts are valid for this block, we free all memory and continue with the next block */
963  if ( ! anycutvalid )
964  {
965  SCIPfreeBufferArray(scip, &ndiagvars);
966  SCIPfreeBufferArray(scip, &diagvars);
967  for (j = blocksize - 1; j >= 0; j--)
968  {
969  SCIPfreeBufferArray(scip, &constnonzeroentries[j]);
970  }
971  SCIPfreeBufferArray(scip, &constnonzeroentries);
972  SCIPfreeBufferArray(scip, &nconstnonzeroentries);
973  continue;
974  }
975 
976  /* find all variables with corresponding diagonal entries for a row with nonzero non-diagonal constant entry, also check for entries
977  * that prevent the cut from being valid */
978  for (v = 0; v < nvars; v++)
979  {
980  for (j = 0; j < consdata->nvarnonz[v]; j++)
981  {
982  /* if it is a diagonal entry for an index that might have a valid cut, we add the variable to the corresponding array if it
983  * is an integer variable and mark the cut invalid otherwise */
984  if ( (consdata->col[v][j] == consdata->row[v][j]) && (nconstnonzeroentries[consdata->col[v][j]] > 0) )
985  {
986  if ( SCIPvarIsIntegral(consdata->vars[v]) )
987  {
988  assert( ! SCIPisEQ(scip, consdata->val[v][j], 0.0) );
989  diagvars[consdata->col[v][j]][ndiagvars[consdata->col[v][j]]] = v;
990  ndiagvars[consdata->col[v][j]]++;
991  }
992  else
993  {
994  nconstnonzeroentries[consdata->col[v][j]] = -2;
995  }
996  }
997  /* If it is a non-diagonal entry, we can no longer use this entry for a cut. If the last entry is removed for a column/row,
998  * mark this column/row invalid (but we still have to free memory later, so we have to set it to -2 instead of 0) */
999  else if ( consdata->col[v][j] != consdata->row[v][j] )
1000  {
1001  if ( nconstnonzeroentries[consdata->col[v][j]] > 0 )
1002  {
1003  /* search for the corresponding row-entry in constnonzeroentries */
1004  for (k = 0; k < nconstnonzeroentries[consdata->col[v][j]]; k++)
1005  {
1006  if ( constnonzeroentries[consdata->col[v][j]][k] == consdata->row[v][j] )
1007  {
1008  /* if there are remaining entries, we shift them back */
1009  if ( nconstnonzeroentries[consdata->col[v][j]] > k + 1 )
1010  {
1011  for (l = k + 1; l < nconstnonzeroentries[consdata->col[v][j]]; l++)
1012  constnonzeroentries[consdata->col[v][j]][l - 1] = constnonzeroentries[consdata->col[v][j]][l];
1013  }
1014  nconstnonzeroentries[consdata->col[v][j]]--;
1015  /* if this was the last entry for this index, we mark it invalid */
1016  if ( nconstnonzeroentries[consdata->col[v][j]] == 0 )
1017  nconstnonzeroentries[consdata->col[v][j]] = -2;
1018  break; /* we should not have another entry for this combination of row and column */
1019  }
1020  }
1021  }
1022  /* do the same for the row */
1023  if ( nconstnonzeroentries[consdata->row[v][j]] > 0 )
1024  {
1025  /* search for the corresponding row-entry in constnonzeroentries */
1026  for (k = 0; k < nconstnonzeroentries[consdata->row[v][j]]; k++)
1027  {
1028  if ( constnonzeroentries[consdata->row[v][j]][k] == consdata->col[v][j] )
1029  {
1030  /* if there are remaining entries, we shift them back */
1031  if ( nconstnonzeroentries[consdata->row[v][j]] > k + 1 )
1032  {
1033  for (l = k + 1; l < nconstnonzeroentries[consdata->row[v][j]]; l++)
1034  constnonzeroentries[consdata->row[v][j]][l - 1] = constnonzeroentries[consdata->row[v][j]][l];
1035  }
1036  nconstnonzeroentries[consdata->row[v][j]]--;
1037  /* if this was the last entry for this index, we mark it invalid */
1038  if ( nconstnonzeroentries[consdata->row[v][j]] == 0 )
1039  nconstnonzeroentries[consdata->row[v][j]] = -2;
1040  break; /* we should not have another entry for this combination of row and column */
1041  }
1042  }
1043  }
1044  }
1045  }
1046  }
1047 
1048  for (j = 0; j < blocksize; ++j)
1049  {
1050  if ( nconstnonzeroentries[j] > 0 )
1051  {
1052  SCIP_CALL( SCIPallocBufferArray(scip, &vals, ndiagvars[j]) );
1053  SCIP_CALL( SCIPallocBufferArray(scip, &vars, ndiagvars[j]) );
1054 
1055  /* get the corresponding SCIP variables and set all coefficients to 1 */
1056  for (v = 0; v < ndiagvars[j]; ++v)
1057  {
1058  vars[v] = consdata->vars[diagvars[j][v]];
1059  vals[v] = 1.0;
1060  }
1061 #ifndef NDEBUG
1062  snprintfreturn = SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "diag_zero_impl_block_%d_row_%d", i, j);
1063  assert( snprintfreturn < SCIP_MAXSTRLEN ); /* check whether name fits into string */
1064 #else
1065  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "diag_zero_impl_block_%d_row_%d", i, j);
1066 #endif
1067 
1068  /* add the linear constraint sum_j 1.0 * diagvars[j] >= 1.0 */
1069  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, cutname, ndiagvars[j], vars, vals, 1.0, SCIPinfinity(scip), TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE) );
1070  SCIP_CALL( SCIPaddCons(scip, cons) );
1071 #ifdef SCIP_MORE_DEBUG
1072  SCIPinfoMessage(scip, NULL, "Added lp-constraint: ");
1073  SCIP_CALL( SCIPprintCons(scip, cons, NULL) );
1074  SCIPinfoMessage(scip, NULL, "\n");
1075 #endif
1076  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
1077  (*naddconss)++;
1078 
1079  SCIPfreeBufferArray(scip, &vars);
1080  SCIPfreeBufferArray(scip, &vals);
1081  }
1082  if ( nconstnonzeroentries[j] == -2 || nconstnonzeroentries[j] > 0 )
1083  {
1084  SCIPfreeBufferArray(scip, &diagvars[j]);
1085  }
1086  }
1087 
1088  SCIPfreeBufferArray(scip, &ndiagvars);
1089  SCIPfreeBufferArray(scip, &diagvars);
1090  for (j = blocksize - 1; j >= 0; j--)
1091  {
1092  SCIPfreeBufferArray(scip, &constnonzeroentries[j]);
1093  }
1094  SCIPfreeBufferArray(scip, &constnonzeroentries);
1095  SCIPfreeBufferArray(scip, &nconstnonzeroentries);
1096  }
1097 
1098  return SCIP_OKAY;
1099 }
1100 
1115 static
1116 SCIP_RETCODE addTwoMinorLinConstraints(
1117  SCIP* scip,
1118  SCIP_CONS** conss,
1119  int nconss,
1120  int* naddconss
1121  )
1122 {
1123  char name[SCIP_MAXSTRLEN];
1124  SCIP_VAR** consvars;
1125  SCIP_Real* consvals;
1126  int blocksize;
1127  int nvars;
1128  int c;
1129  int i;
1130 
1131  assert( scip != NULL );
1132  assert( naddconss != NULL );
1133 
1134  for (c = 0; c < nconss; ++c)
1135  {
1136  SCIP_CONSDATA* consdata;
1137  SCIP_Real** matrices = NULL;
1138  SCIP_Real* constmatrix;
1139  SCIP_Real* coef;
1140  int s;
1141  int t;
1142 
1143  assert( conss[c] != NULL );
1144  consdata = SCIPconsGetData(conss[c]);
1145  assert( consdata != NULL );
1146 
1147  blocksize = consdata->blocksize;
1148  nvars = consdata->nvars;
1149 
1150  SCIP_CALL( SCIPallocBufferArray(scip, &consvars, nvars) );
1151  SCIP_CALL( SCIPallocBufferArray(scip, &consvals, nvars) );
1152  SCIP_CALL( SCIPallocBufferArray(scip, &coef, nvars) );
1153 
1154  /* get matrices */
1155  SCIP_CALL( SCIPallocBufferArray(scip, &constmatrix, blocksize * blocksize) );
1156  SCIP_CALL( SCIPconsSdpGetFullConstMatrix(scip, conss[c], constmatrix) );
1157 
1158  SCIP_CALL( SCIPallocBufferArray(scip, &matrices, nvars) );
1159  for (i = 0; i < nvars; ++i)
1160  {
1161  SCIP_CALL( SCIPallocBufferArray(scip, &matrices[i], blocksize * blocksize) );
1162  SCIP_CALL( SCIPconsSdpGetFullAj(scip, conss[c], i, matrices[i]) );
1163  }
1164 
1165  /* loop over all possible entries */
1166  for (s = 0; s < blocksize; ++s)
1167  {
1168  for (t = 0; t < s; ++t)
1169  {
1170  SCIP_CONS* cons;
1171  SCIP_Real val;
1172  SCIP_Real lhs;
1173  SCIP_Real activitylb = 0.0;
1174  int nconsvars = 0;
1175  int cnt = 0;
1176 
1177  /* skip diagonal entries */
1178  if ( s == t )
1179  continue;
1180 
1181  /* collect coefficients */
1182  BMSclearMemoryArray(coef, nvars);
1183  for (i = 0; i < nvars; ++i)
1184  {
1185  val = matrices[i][s * blocksize + t];
1186  if ( ! SCIPisZero(scip, val) )
1187  {
1188  coef[i] = -2.0 * val;
1189  ++cnt;
1190  }
1191  }
1192 
1193  /* only proceed if off-diagonal is nonempty */
1194  if ( cnt == 0 )
1195  continue;
1196 
1197  /* add diagonal entries for s */
1198  for (i = 0; i < nvars; ++i)
1199  {
1200  val = matrices[i][s * blocksize + s];
1201  if ( ! SCIPisZero(scip, val) )
1202  coef[i] += val;
1203  }
1204 
1205  /* add diagonal entries for t */
1206  for (i = 0; i < nvars; ++i)
1207  {
1208  val = matrices[i][t * blocksize + t];
1209  if ( ! SCIPisZero(scip, val) )
1210  coef[i] += val;
1211  }
1212 
1213  /* get constraint */
1214  for (i = 0; i < nvars; ++i)
1215  {
1216  if ( ! SCIPisZero(scip, coef[i]) )
1217  {
1218  consvals[nconsvars] = coef[i];
1219  consvars[nconsvars] = consdata->vars[i];
1220  ++nconsvars;
1221 
1222  /* compute lower bound on activity */
1223  if ( coef[i] > 0 )
1224  activitylb += coef[i] * SCIPvarGetLbGlobal(consdata->vars[i]);
1225  else
1226  activitylb += coef[i] * SCIPvarGetUbGlobal(consdata->vars[i]);
1227  }
1228  }
1229 
1230  /* only proceed if cut is nontrivial */
1231  if ( nconsvars <= 0 )
1232  continue;
1233 
1234  /* compute rhs */
1235  lhs = constmatrix[s * blocksize + s] + constmatrix[t * blocksize + t] - 2.0 * constmatrix[s * blocksize + t];
1236 
1237  /* only proceed if constraint is not redundant */
1238  if ( SCIPisGE(scip, activitylb, lhs) )
1239  continue;
1240 
1241  /* add linear constraint (only propagate) */
1242  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "2x2minorlin#%d#%d", s, t);
1243  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, name, nconsvars, consvars, consvals, lhs, SCIPinfinity(scip),
1244  TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE) );
1245  SCIP_CALL( SCIPaddCons(scip, cons) );
1246 #ifdef SCIP_MORE_DEBUG
1247  SCIPinfoMessage(scip, NULL, "Added 2x2 minor linear constraint: ");
1248  SCIP_CALL( SCIPprintCons(scip, cons, NULL) );
1249  SCIPinfoMessage(scip, NULL, "\n");
1250 #endif
1251  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
1252  ++(*naddconss);
1253  }
1254  }
1255 
1256  for (i = 0; i < nvars; ++i)
1257  SCIPfreeBufferArray(scip, &matrices[i]);
1258  SCIPfreeBufferArray(scip, &matrices);
1259  SCIPfreeBufferArray(scip, &constmatrix);
1260  SCIPfreeBufferArray(scip, &coef);
1261  SCIPfreeBufferArray(scip, &consvals);
1262  SCIPfreeBufferArray(scip, &consvars);
1263  }
1264 
1265  return SCIP_OKAY;
1266 }
1267 
1273 static
1274 SCIP_RETCODE addTwoMinorProdConstraints(
1275  SCIP* scip,
1276  SCIP_CONS** conss,
1277  int nconss,
1278  int* naddconss
1279  )
1280 {
1281  char name[SCIP_MAXSTRLEN];
1282  SCIP_VAR** consvars;
1283  SCIP_Real* consvals;
1284  int blocksize;
1285  int nvars;
1286  int c;
1287  int i;
1288  int j;
1289 
1290  assert( scip != NULL );
1291  assert( naddconss != NULL );
1292 
1293  for (c = 0; c < nconss; ++c)
1294  {
1295  SCIP_CONSDATA* consdata;
1296  SCIP_Real** matrices = NULL;
1297  SCIP_Real* constmatrix;
1298 
1299  assert( conss[c] != NULL );
1300  consdata = SCIPconsGetData(conss[c]);
1301  assert( consdata != NULL );
1302 
1303  blocksize = consdata->blocksize;
1304  nvars = consdata->nvars;
1305 
1306  SCIP_CALL( SCIPallocBufferArray(scip, &consvars, nvars) );
1307  SCIP_CALL( SCIPallocBufferArray(scip, &consvals, nvars) );
1308 
1309  SCIP_CALL( SCIPallocBufferArray(scip, &constmatrix, blocksize * blocksize) );
1310  SCIP_CALL( SCIPconsSdpGetFullConstMatrix(scip, conss[c], constmatrix) );
1311 
1312  /* loop over all nonzero matrix entries in the constant part */
1313  for (j = 0; j < consdata->constnnonz; ++j)
1314  {
1315  SCIP_Real val;
1316  SCIP_Real prod;
1317  int s;
1318  int t;
1319 
1320  s = consdata->constrow[j];
1321  t = consdata->constcol[j];
1322 
1323  /* skip diagonal entries */
1324  if ( s == t )
1325  continue;
1326 
1327  /* compute matrices if not yet done */
1328  if ( matrices == NULL )
1329  {
1330  SCIP_CALL( SCIPallocBufferArray(scip, &matrices, nvars) );
1331  for (i = 0; i < nvars; ++i)
1332  {
1333  SCIP_CALL( SCIPallocBufferArray(scip, &matrices[i], blocksize * blocksize) );
1334  SCIP_CALL( SCIPconsSdpGetFullAj(scip, conss[c], i, matrices[i]) );
1335  }
1336  }
1337  assert( matrices != NULL );
1338 
1339  /* check whether diagonal entries in the matrices are all 0 */
1340  for (i = 0; i < nvars; ++i)
1341  {
1342  if ( ! SCIPisZero(scip, matrices[i][s * blocksize + s]) )
1343  break;
1344  if ( ! SCIPisZero(scip, matrices[i][t * blocksize + t]) )
1345  break;
1346  }
1347  if ( i < nvars )
1348  continue;
1349 
1350  /* check if product is sufficiently positive */
1351  prod = constmatrix[s * blocksize + s] * constmatrix[t * blocksize + t];
1352  if ( SCIPisEfficacious(scip, prod) )
1353  {
1354  SCIP_CONS* cons;
1355  SCIP_Real activitylb = 0.0;
1356  SCIP_Real lhs;
1357  int nconsvars = 0;
1358 
1359  /* fill in constraint */
1360  for (i = 0; i < nvars; ++i)
1361  {
1362  val = matrices[i][s * blocksize + t];
1363  if ( ! SCIPisZero(scip, val) )
1364  {
1365  consvals[nconsvars] = val;
1366  consvars[nconsvars] = consdata->vars[i];
1367  ++nconsvars;
1368 
1369  /* compute lower bound on activity */
1370  if ( val > 0 )
1371  activitylb += val * SCIPvarGetLbGlobal(consdata->vars[i]);
1372  else
1373  activitylb += val * SCIPvarGetUbGlobal(consdata->vars[i]);
1374  }
1375  }
1376 
1377  lhs = consdata->constval[j] - sqrt(prod);
1378 
1379  /* only proceed if constraint is not redundant */
1380  if ( SCIPisGE(scip, activitylb, lhs) )
1381  continue;
1382 
1383  /* add linear constraint (only propagate) */
1384  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "2x2minorprod#%d#%d", s, t);
1385  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, name, nconsvars, consvars, consvals, lhs, SCIPinfinity(scip),
1386  TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE) );
1387  SCIP_CALL( SCIPaddCons(scip, cons) );
1388 #ifdef SCIP_MORE_DEBUG
1389  SCIPinfoMessage(scip, NULL, "Added 2x2 minor product constraint: ");
1390  SCIP_CALL( SCIPprintCons(scip, cons, NULL) );
1391  SCIPinfoMessage(scip, NULL, "\n");
1392 #endif
1393  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
1394  ++(*naddconss);
1395  }
1396  }
1397 
1398  if ( matrices != NULL )
1399  {
1400  for (i = 0; i < nvars; ++i)
1401  SCIPfreeBufferArray(scip, &matrices[i]);
1402  SCIPfreeBufferArray(scip, &matrices);
1403  }
1404  SCIPfreeBufferArray(scip, &constmatrix);
1405  SCIPfreeBufferArray(scip, &consvals);
1406  SCIPfreeBufferArray(scip, &consvars);
1407  }
1408 
1409  return SCIP_OKAY;
1410 }
1411 
1413 static
1414 SCIP_RETCODE move_1x1_blocks_to_lp(
1415  SCIP* scip,
1416  SCIP_CONSHDLR* conshdlr,
1417  SCIP_CONS** conss,
1418  int nconss,
1419  int* naddconss,
1420  int* ndelconss,
1421  int* nchgbds,
1422  SCIP_RESULT* result
1423  )
1424 {
1425  char cutname[SCIP_MAXSTRLEN];
1426  SCIP_CONSDATA* consdata;
1427  SCIP_CONSHDLRDATA* conshdlrdata;
1428  SCIP_CONS* cons;
1429  SCIP_VAR** vars;
1430  SCIP_Real* coeffs;
1431  SCIP_Real rhs;
1432  int nnonz;
1433  int nvars;
1434  int i;
1435  int j;
1436  int v;
1437 #ifndef NDEBUG
1438  int snprintfreturn; /* used to assert the return code of snprintf */
1439 #endif
1440 
1441  assert( scip != NULL );
1442  assert( conshdlr != NULL );
1443  assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 || strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLRRANK1_NAME) == 0 );
1444  assert( result != NULL );
1445 
1446  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1447  assert( conshdlrdata != NULL );
1448 
1449  for (i = 0; i < nconss && *result != SCIP_CUTOFF; ++i)
1450  {
1451  consdata = SCIPconsGetData(conss[i]);
1452  assert( consdata != NULL );
1453 
1454  /* if there is a 1x1 SDP-Block */
1455  if ( consdata->blocksize == 1 )
1456  {
1457  int cnt = 0;
1458 
1459  nvars = consdata->nvars;
1460  nnonz = consdata->nnonz;
1461  SCIP_CALL( SCIPallocBufferArray(scip, &vars, nvars) );
1462  SCIP_CALL( SCIPallocBufferArray(scip, &coeffs, nnonz) );
1463 
1464  /* get all lhs-entries */
1465  for (v = 0; v < nvars; v++)
1466  {
1467  assert( consdata->nvarnonz[v] <= 1 ); /* in a 1x1 block there may be at most one entry per variable */
1468 
1469  for (j = 0; j < consdata->nvarnonz[v]; j++)
1470  {
1471  assert( consdata->col[v][j] == 0 && consdata->row[v][j] == 0 ); /* if the block has size 1, all entries should have row and col equal to 0 */
1472  if ( ! SCIPisZero(scip, consdata->val[v][j]) )
1473  {
1474  coeffs[cnt] = consdata->val[v][j];
1475  vars[cnt] = consdata->vars[v];
1476  ++cnt;
1477  }
1478  }
1479  }
1480 
1481  /* get rhs */
1482  assert( consdata->constnnonz <= 1 ); /* the 1x1 constant matrix may only have one entry */
1483 
1484  rhs = (consdata->constnnonz == 1) ? consdata->constval[0] : 0.0; /* if this one entry is not 0, than this is the rhs, otherwise it's 0 */
1485 
1486  /* if there is more than one nonzero left-hand-side-entry, add a linear constraint, otherwise update the variable bound */
1487  if ( cnt > 1 )
1488  {
1489  /* add new linear cons */
1490 #ifndef NDEBUG
1491  snprintfreturn = SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "1x1block_%d", ++(conshdlrdata->n1x1blocks));
1492  assert( snprintfreturn < SCIP_MAXSTRLEN ); /* check whether name fits into string */
1493 #else
1494  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "1x1block_%d", ++(conshdlrdata->n1x1blocks));
1495 #endif
1496 
1497  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, cutname, cnt, vars, coeffs, rhs, SCIPinfinity(scip),
1498  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE) );
1499 
1500  SCIP_CALL( SCIPaddCons(scip, cons) );
1501 #ifdef SCIP_MORE_DEBUG
1502  SCIPinfoMessage(scip, NULL, "Added lp-constraint: ");
1503  SCIP_CALL( SCIPprintCons(scip, cons, NULL) );
1504  SCIPinfoMessage(scip, NULL, "\n");
1505 #endif
1506  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
1507 
1508  (*naddconss)++;
1509  }
1510  else if ( cnt == 1 )
1511  {
1512  SCIP_Bool infeasible = FALSE;
1513  SCIP_Bool tightened;
1514 
1515  /* try to tighten bound */
1516  if ( SCIPisPositive(scip, coeffs[0]) )
1517  {
1518  SCIP_CALL( SCIPtightenVarLb(scip, vars[0], rhs / coeffs[0], FALSE, &infeasible, &tightened) );
1519  if ( tightened )
1520  {
1521  SCIPdebugMsg(scip, "Tightend lower bound of <%s> to %g because of diagonal values of SDP-constraint %s!\n",
1522  SCIPvarGetName(vars[0]), SCIPvarGetLbGlobal(vars[0]), SCIPconsGetName(conss[i]));
1523  ++(*nchgbds);
1524  }
1525  }
1526  else
1527  {
1528  assert( SCIPisNegative(scip, coeffs[0]) );
1529  SCIP_CALL( SCIPtightenVarUb(scip, vars[0], rhs / coeffs[0], FALSE, &infeasible, &tightened) );
1530  if ( tightened )
1531  {
1532  SCIPdebugMsg(scip, "Tightend upper bound of <%s> to %g because of diagonal values of SDP-constraint %s!\n",
1533  SCIPvarGetName(vars[0]), SCIPvarGetUbGlobal(vars[0]), SCIPconsGetName(conss[i]));
1534  ++(*nchgbds);
1535  }
1536  }
1537 
1538  if ( infeasible )
1539  *result = SCIP_CUTOFF;
1540  }
1541  else
1542  {
1543  assert( cnt == 0 );
1544  SCIPdebugMsg(scip, "Detected 1x1 SDP-block without any nonzero coefficients \n");
1545  if ( SCIPisFeasGT(scip, rhs, 0.0) )
1546  {
1547  SCIPdebugMsg(scip, "Detected infeasibility in 1x1 SDP-block without any nonzero coefficients but with strictly positive rhs\n");
1548  *result = SCIP_CUTOFF;
1549  }
1550  }
1551 
1552  /* delete old 1x1 sdpcone */
1553  SCIP_CALL( SCIPdelCons(scip, conss[i]) );
1554  (*ndelconss)++;
1555 
1556  SCIPfreeBufferArray(scip, &coeffs);
1557  SCIPfreeBufferArray(scip, &vars);
1558  }
1559  }
1560  return SCIP_OKAY;
1561 }
1562 
1564 static
1565 SCIP_RETCODE unlockVar(
1566  SCIP* scip,
1567  SCIP_CONSDATA* consdata,
1568  int v
1569  )
1570 {
1571  assert( scip != NULL );
1572  assert( consdata != NULL );
1573  assert( 0 <= v && v < consdata->nvars );
1574 
1575  if ( consdata->locks != NULL )
1576  {
1577  assert( consdata->locks[v] == -2 || consdata->locks[v] == -1 || consdata->locks[v] == 0 || consdata->locks[v] == 1 );
1578 
1579  if ( consdata->locks[v] == 1 )
1580  {
1581  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], SCIP_LOCKTYPE_MODEL, 0, -1) );
1582  }
1583  else if ( consdata->locks[v] == - 1 )
1584  {
1585  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], SCIP_LOCKTYPE_MODEL, -1, 0) );
1586  }
1587  else if ( consdata->locks[v] == 0 )
1588  {
1589  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], SCIP_LOCKTYPE_MODEL, -1, -1) );
1590  }
1591  }
1592 
1593  return SCIP_OKAY;
1594 }
1595 
1597 static
1598 SCIP_RETCODE updateVarLocks(
1599  SCIP* scip,
1600  SCIP_CONS* cons,
1601  int v
1602  )
1603 {
1604  SCIP_CONSDATA* consdata;
1605  SCIP_Real eigenvalue;
1606  SCIP_Real* Aj;
1607  int blocksize;
1608  int newlock = -2;
1609 
1610  assert( scip != NULL );
1611  assert( cons != NULL );
1612 
1613  consdata = SCIPconsGetData(cons);
1614  assert( consdata != NULL );
1615  assert( consdata->locks != NULL );
1616  assert( 0 <= v && v < consdata->nvars );
1617 
1618  /* rank-1 constraints are always up- and down-locked */
1619  if ( consdata->rankone )
1620  {
1621  SCIP_CALL( unlockVar(scip, consdata, v) );
1622  consdata->locks[v] = 0;
1623  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], SCIP_LOCKTYPE_MODEL, 1, 1) );
1624  return SCIP_OKAY;
1625  }
1626 
1627  blocksize = consdata->blocksize;
1628  SCIP_CALL( SCIPallocBufferArray(scip, &Aj, blocksize * blocksize) );
1629 
1630  SCIP_CALL( SCIPconsSdpGetFullAj(scip, cons, v, Aj) );
1631 
1632  /* compute new lock as in consLockSdp */
1633  SCIP_CALL( SCIPlapackComputeIthEigenvalue(SCIPbuffer(scip), FALSE, blocksize, Aj, 1, &eigenvalue, NULL) );
1634  if ( SCIPisNegative(scip, eigenvalue) )
1635  newlock = 1; /* up-lock */
1636 
1637  if ( SCIPisPositive(scip, eigenvalue) )
1638  newlock = -1; /* down-lock */
1639  else
1640  {
1641  SCIP_CALL( SCIPlapackComputeIthEigenvalue(SCIPbuffer(scip), FALSE, blocksize, Aj, blocksize, &eigenvalue, NULL) );
1642  if ( SCIPisPositive(scip, eigenvalue) )
1643  {
1644  if ( newlock == 1 )
1645  newlock = 0; /* up- and down-lock */
1646  else
1647  newlock = -1; /* down-lock */
1648  }
1649  }
1650 
1651  SCIPfreeBufferArray(scip, &Aj);
1652 
1653  /* if new lock is not equal to the old one, unlock variable and add new locks */
1654  if ( newlock != consdata->locks[v] )
1655  {
1656  SCIP_CALL( unlockVar(scip, consdata, v) );
1657  if ( newlock == 1 ) /* up-lock */
1658  {
1659  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], SCIP_LOCKTYPE_MODEL, 0, 1) );
1660  }
1661  else if ( newlock == -1 ) /* down-lock */
1662  {
1663  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], SCIP_LOCKTYPE_MODEL, 1, 0) );
1664  }
1665  else if ( newlock == 0 ) /* up and down lock */
1666  {
1667  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], SCIP_LOCKTYPE_MODEL, 1, 1) );
1668  }
1669  else
1670  assert( newlock == -2 );
1671 
1672  consdata->locks[v] = newlock;
1673  }
1674 
1675  return SCIP_OKAY;
1676 }
1677 
1678 #ifndef NDEBUG
1679 
1680 static
1681 SCIP_RETCODE checkVarsLocks(
1682  SCIP* scip,
1683  SCIP_CONS* cons
1684  )
1685 {
1686  SCIP_CONSDATA* consdata;
1687  SCIP_Real* Aj;
1688  int blocksize;
1689  int v;
1690 
1691  assert( scip != NULL );
1692  assert( cons != NULL );
1693 
1694  consdata = SCIPconsGetData(cons);
1695  assert( consdata != NULL );
1696  assert( consdata->locks != NULL );
1697 
1698  /* rank-1 constraints should always be up- and down-locked */
1699  if ( consdata->rankone )
1700  {
1701  for (v = 0; v < consdata->nvars; ++v)
1702  assert( consdata->locks[v] == 0 );
1703  return SCIP_OKAY;
1704  }
1705 
1706  blocksize = consdata->blocksize;
1707  SCIP_CALL( SCIPallocBufferArray(scip, &Aj, blocksize * blocksize) );
1708 
1709  for (v = 0; v < consdata->nvars; ++v)
1710  {
1711  SCIP_Real eigenvalue;
1712  int newlock = -2;
1713 
1714  SCIP_CALL( SCIPconsSdpGetFullAj(scip, cons, v, Aj) );
1715 
1716  /* compute new lock as in consLockSdp */
1717  SCIP_CALL( SCIPlapackComputeIthEigenvalue(SCIPbuffer(scip), FALSE, blocksize, Aj, 1, &eigenvalue, NULL) );
1718  if ( SCIPisNegative(scip, eigenvalue) )
1719  newlock = 1; /* up-lock */
1720 
1721  if ( SCIPisPositive(scip, eigenvalue) )
1722  newlock = -1; /* down-lock */
1723  else
1724  {
1725  SCIP_CALL( SCIPlapackComputeIthEigenvalue(SCIPbuffer(scip), FALSE, blocksize, Aj, blocksize, &eigenvalue, NULL) );
1726  if ( SCIPisPositive(scip, eigenvalue) )
1727  {
1728  if ( newlock == 1 )
1729  newlock = 0; /* up- and down-lock */
1730  else
1731  newlock = -1; /* down-lock */
1732  }
1733  }
1734 
1735  assert( newlock == consdata->locks[v] );
1736  }
1737 
1738  SCIPfreeBufferArray(scip, &Aj);
1739 
1740  return SCIP_OKAY;
1741 }
1742 #endif
1743 
1745 static
1746 SCIP_RETCODE multiaggrVar(
1747  SCIP* scip,
1748  SCIP_CONS* cons,
1749  int v,
1750  SCIP_VAR** aggrvars,
1751  SCIP_Real* scalars,
1752  int naggrvars,
1753  SCIP_Real constant,
1754  int* savedcol,
1755  int* savedrow,
1756  SCIP_Real* savedval,
1757  int* nfixednonz,
1758  int* vararraylength
1759  )
1760 {
1761  SCIP_CONSDATA* consdata;
1762  int startind;
1763  int aggrind;
1764  int aggrtargetlength;
1765  int globalnvars;
1766  int aggrconsind;
1767  int i;
1768 
1769  assert( scip != NULL );
1770  assert( cons != NULL );
1771  assert( scalars != NULL );
1772  assert( naggrvars > 0 );
1773  assert( savedcol != NULL );
1774  assert( savedrow != NULL );
1775  assert( savedval != NULL );
1776  assert( nfixednonz != NULL );
1777 
1778  consdata = SCIPconsGetData(cons);
1779  assert( consdata != NULL );
1780  assert( consdata->locks != NULL );
1781  assert( 0 <= v && v < consdata->nvars );
1782 
1783  /* save the current nfixednonz-index, all entries starting from here will need to be added to the variables this is aggregated to */
1784  startind = *nfixednonz;
1785 
1786  if ( SCIPisEQ(scip, constant, 0.0) )
1787  {
1788  /* if there is no constant part, we just save the nonzeros to add them to the variables they are aggregated to, we do this to be able to remove
1789  * the nonzero-arrays for this variable to be able to fill it with a newly inserted variable, as copying all variables, if we created an empty
1790  * gap in the variable arrays, will be more time consuming then copying all variables (as we will usually have more variables than nonzeros per
1791  * variable */
1792  for (i = 0; i < consdata->nvarnonz[v]; i++)
1793  {
1794  savedcol[*nfixednonz] = consdata->col[v][i];
1795  savedrow[*nfixednonz] = consdata->row[v][i];
1796  savedval[*nfixednonz] = consdata->val[v][i];
1797  (*nfixednonz)++;
1798  }
1799  }
1800  else
1801  {
1802  for (i = 0; i < consdata->nvarnonz[v]; i++)
1803  {
1804  savedcol[*nfixednonz] = consdata->col[v][i];
1805  savedrow[*nfixednonz] = consdata->row[v][i];
1806  savedval[*nfixednonz] = consdata->val[v][i] * constant;
1807  (*nfixednonz)++;
1808  }
1809  }
1810  assert( *nfixednonz - startind == consdata->nvarnonz[v] );
1811 
1812  /* sort them by nondecreasing row and then col to make the search for already existing entries easier (this is done here, because it
1813  * only needs to be done once and not for each variable this is multiaggregated to), we add startind to the pointers to only start where we started
1814  * inserting, the number of elements added to the saved arrays for this variable is nfixednonz - startind */
1815  SCIPsdpVarfixerSortRowCol(savedrow + startind, savedcol + startind, savedval + startind, *nfixednonz - startind);
1816 
1817  /* free the memory for the entries of the aggregated variable */
1818  SCIPfreeBlockMemoryArray(scip, &(consdata->val[v]), consdata->nvarnonz[v]);
1819  SCIPfreeBlockMemoryArray(scip, &(consdata->row[v]), consdata->nvarnonz[v]);
1820  SCIPfreeBlockMemoryArray(scip, &(consdata->col[v]), consdata->nvarnonz[v]);
1821 
1822  /* unlock variable */
1823  SCIP_CALL( unlockVar(scip, consdata, v) );
1824 
1825  /* fill the empty spot of the (multi-)aggregated variable with the last variable of this constraint (as they don't have to be sorted) */
1826  SCIP_CALL( SCIPreleaseVar(scip, &consdata->vars[v]) );
1827  consdata->col[v] = consdata->col[consdata->nvars - 1];
1828  consdata->row[v] = consdata->row[consdata->nvars - 1];
1829  consdata->val[v] = consdata->val[consdata->nvars - 1];
1830  consdata->nvarnonz[v] = consdata->nvarnonz[consdata->nvars - 1];
1831  consdata->vars[v] = consdata->vars[consdata->nvars - 1];
1832  consdata->locks[v] = consdata->locks[consdata->nvars - 1];
1833  (consdata->nvars)--;
1834 
1835  /* iterate over all variables this was aggregated to and insert the corresponding nonzeros */
1836  for (aggrind = 0; aggrind < naggrvars; aggrind++)
1837  {
1838  /* check if the variable already exists in this block */
1839  aggrconsind = -1;
1840  for (i = 0; i < consdata->nvars; i++)
1841  {
1842  if ( consdata->vars[i] == aggrvars[aggrind] )
1843  {
1844  aggrconsind = i;
1845  break;
1846  }
1847  }
1848 
1849  if ( aggrconsind > -1 )
1850  {
1851  /* if the varialbe to aggregate to is already part of this sdp-constraint, just add the nonzeros of the old variable to it */
1852 
1853  /* resize the arrays to the maximally needed length */
1854  aggrtargetlength = consdata->nvarnonz[aggrconsind] + *nfixednonz - startind;
1855  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(consdata->row[aggrconsind]), consdata->nvarnonz[aggrconsind], aggrtargetlength) );
1856  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(consdata->col[aggrconsind]), consdata->nvarnonz[aggrconsind], aggrtargetlength) );
1857  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(consdata->val[aggrconsind]), consdata->nvarnonz[aggrconsind], aggrtargetlength) );
1858 
1859  if ( SCIPisEQ(scip, constant, 0.0) )
1860  {
1861  /* in this case we saved the original values in savedval, we add startind to the pointers to only add those from
1862  * the current variable, the number of entries is the current position minus the position whre we started */
1863  SCIP_CALL( SCIPsdpVarfixerMergeArrays(SCIPblkmem(scip), SCIPepsilon(scip), savedrow + startind, savedcol + startind, savedval + startind,
1864  *nfixednonz - startind, TRUE, scalars[aggrind], consdata->row[aggrconsind], consdata->col[aggrconsind],
1865  consdata->val[aggrconsind], &(consdata->nvarnonz[aggrconsind]), aggrtargetlength) );
1866  }
1867  else
1868  {
1869  /* in this case we saved the original values * constant, so we now have to divide by constant, we add startind to the pointers
1870  * to only add those from the current variable, the number of entries is the current position minus the position whre we started */
1871  SCIP_CALL( SCIPsdpVarfixerMergeArrays(SCIPblkmem(scip), SCIPepsilon(scip), savedrow + startind, savedcol + startind, savedval + startind,
1872  *nfixednonz - startind, TRUE, scalars[aggrind] / constant, consdata->row[aggrconsind], consdata->col[aggrconsind],
1873  consdata->val[aggrconsind], &(consdata->nvarnonz[aggrconsind]), aggrtargetlength) );
1874  }
1875 
1876  /* shrink them again if nonzeros could be combined */
1877  assert( consdata->nvarnonz[aggrconsind] <= aggrtargetlength );
1878  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(consdata->row[aggrconsind]), aggrtargetlength, consdata->nvarnonz[aggrconsind]) );
1879  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(consdata->col[aggrconsind]), aggrtargetlength, consdata->nvarnonz[aggrconsind]) );
1880  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(consdata->val[aggrconsind]), aggrtargetlength, consdata->nvarnonz[aggrconsind]) );
1881 
1882  SCIP_CALL( updateVarLocks(scip, cons, aggrconsind) );
1883  }
1884  else
1885  {
1886  /* the variable has to be added to this constraint */
1887  SCIPdebugMsg(scip, "adding variable %s to SDP constraint %s because of (multi-)aggregation\n", SCIPvarGetName(aggrvars[aggrind]), SCIPconsGetName(cons));
1888 
1889  /* check if we have to enlarge the arrays */
1890  if ( consdata->nvars == *vararraylength )
1891  {
1892  globalnvars = SCIPgetNVars(scip);
1893 
1894  /* we don't want to enlarge this by one for every variable added, so we immediately set it to the maximum possible size */
1895  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->col, *vararraylength, globalnvars) );
1896  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->row, *vararraylength, globalnvars) );
1897  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->val, *vararraylength, globalnvars) );
1898  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->nvarnonz, *vararraylength, globalnvars) );
1899  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->vars, *vararraylength, globalnvars) );
1900  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->locks, *vararraylength, globalnvars) );
1901  *vararraylength = globalnvars;
1902  }
1903 
1904  /* we insert this variable at the last position, as the ordering doesn't matter */
1905  SCIP_CALL( SCIPcaptureVar(scip, aggrvars[aggrind]) );
1906  consdata->vars[consdata->nvars] = aggrvars[aggrind];
1907 
1908  /* as there were no nonzeros thus far, we can just duplicate the saved arrays to get the nonzeros for the new variable */
1909  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(consdata->col[consdata->nvars]), savedcol + startind, *nfixednonz - startind) ); /*lint !e776*/
1910  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(consdata->row[consdata->nvars]), savedrow + startind, *nfixednonz - startind) ); /*lint !e776*/
1911 
1912  /* if scalars[aggrind] = constant, we would multiply with 1, if constant = 0, we didn't divide by constant, so in these cases, we can just
1913  * memcopy the array of nonzero-values */
1914  /* TODO: only checking scalar and constant for feas-equality might lead to big differences, if the nonzeros they are multiplied with are big,
1915  * e.g. scalar = 0, constant = 10^(-6), nonzero = 10^(10) leads to new nonzero of 10^4 instead of 0 */
1916  if ( SCIPisEQ(scip, scalars[aggrind], constant) || SCIPisEQ(scip, constant, 0.0) )
1917  {
1918  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(consdata->val[consdata->nvars]), savedval + startind, *nfixednonz - startind) ); /*lint !e776*/
1919  consdata->nvarnonz[consdata->nvars] = *nfixednonz - startind; /* as there were no nonzeros thus far, the number of nonzeros equals
1920  * the number of nonzeros of the aggregated variable */
1921  }
1922  else /* we have to multiply all entries by scalar before inserting them */
1923  {
1924  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(consdata->val[consdata->nvars]), *nfixednonz - startind) );
1925 
1926  consdata->nvarnonz[consdata->nvars] = 0;
1927 
1928  for (i = 0; i < *nfixednonz - startind; i++)
1929  {
1930  /* if both scalar and savedval are small this might become too small */
1931  if ( SCIPisPositive(scip, (scalars[i] / constant) * savedval[startind + i]) )
1932  {
1933  consdata->val[consdata->nvars][consdata->nvarnonz[consdata->nvars]] = (scalars[i] / constant) * savedval[startind + i]; /*lint !e679*/
1934  consdata->nvarnonz[consdata->nvars]++;
1935  }
1936  }
1937  }
1938 
1939  consdata->locks[consdata->nvars] = -2;
1940  if ( consdata->nvarnonz[consdata->nvars] > 0 ) /* if scalar and all savedvals were to small */
1941  {
1942  consdata->nvars++;
1943  SCIP_CALL( updateVarLocks(scip, cons, consdata->nvars-1) );
1944  }
1945  else
1946  SCIP_CALL( updateVarLocks(scip, cons, consdata->nvars) );
1947  }
1948  }
1949 
1950  /* if the constant part is zero, we may delete the nonzero-entries from the saved arrays (by resetting the nfixednonz entry to where
1951  * it started, so that these entries will be overwritten */
1952  if ( SCIPisEQ(scip, constant, 0.0) )
1953  *nfixednonz = startind;
1954 
1955 #ifndef NDEBUG
1956  SCIP_CALL( checkVarsLocks(scip, cons) );
1957 #endif
1958 
1959  return SCIP_OKAY;
1960 }
1961 
1962 
1964 static
1965 SCIP_RETCODE fixAndAggrVars(
1966  SCIP* scip,
1967  SCIP_CONS** conss,
1968  int nconss,
1969  SCIP_Bool aggregate
1970  )
1971 {
1972  SCIP_CONSDATA* consdata;
1973  int i;
1974  int* savedcol;
1975  int* savedrow;
1976  SCIP_Real* savedval;
1977  int c;
1978  int v;
1979  int arraylength;
1980  SCIP_VAR* var;
1981  SCIP_VAR** aggrvars;
1982  SCIP_Real scalar;
1983  SCIP_Real* scalars;
1984  int naggrvars;
1985  SCIP_Real constant;
1986  int requiredsize;
1987  int globalnvars;
1988  int vararraylength;
1989 
1990  /* Loop over all variables once, add all fixed to savedrow/col/val; for all multiaggregated variables, if constant-scalar != 0, add
1991  * constant-scalar * entry to savedrow/col/val and call mergeArrays for all aggrvars for savedrow[startindex of this var] and scalar/constant-scalar;
1992  * if constant-scalar == 0, add 1*entry to savedrow/col/val, call mergeArrays for all aggrvars for savedrow[startindex of this var] and scalar and later
1993  * reduce the saved size of savedrow/col/val by the number of nonzeros of the multiaggregated variable to not add them to the constant part later. */
1994 
1995  assert( scip != NULL );
1996  assert( conss != NULL );
1997  assert( nconss >= 0 );
1998 
1999  SCIPdebugMsg(scip, "Calling fixAndAggrVars with aggregate = %u.\n", aggregate);
2000 
2001  for (c = 0; c < nconss; ++c)
2002  {
2003  int nfixednonz = 0;
2004 
2005  assert( conss[c] != NULL );
2006 
2007  consdata = SCIPconsGetData(conss[c]);
2008  assert( consdata != NULL );
2009  assert( consdata->locks != NULL );
2010  assert( consdata->rankone || strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(conss[c])), CONSHDLR_NAME) == 0 );
2011  assert( ! consdata->rankone || strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(conss[c])), CONSHDLRRANK1_NAME) == 0 );
2012 
2013  /* allocate memory to save nonzeros that need to be fixed */
2014  SCIP_CALL( SCIPallocBufferArray(scip, &savedcol, consdata->nnonz) );
2015  SCIP_CALL( SCIPallocBufferArray(scip, &savedrow, consdata->nnonz) );
2016  SCIP_CALL( SCIPallocBufferArray(scip, &savedval, consdata->nnonz) );
2017 
2018  vararraylength = consdata->nvars;
2019  globalnvars = SCIPgetNVars(scip);
2020 
2021  for (v = 0; v < consdata->nvars; v++)/*lint --e{850}*/
2022  {
2023  SCIP_Bool negated = FALSE;
2024 
2025  /* if the variable is negated, get the negation var */
2026  if ( SCIPvarIsBinary(consdata->vars[v]) && SCIPvarIsNegated(consdata->vars[v]) )
2027  {
2028  negated = TRUE;
2029  var = SCIPvarGetNegationVar(consdata->vars[v]);
2030  }
2031  else
2032  var = consdata->vars[v];
2033 
2034  /* check if the variable is fixed in SCIP */
2035  if ( SCIPvarGetStatus(var) == SCIP_VARSTATUS_FIXED || SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var)) )
2036  {
2037  assert( SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var)) );
2038 
2039  SCIPdebugMsg(scip, "Treating globally fixed variable %s with value %f!\n", SCIPvarGetName(var), SCIPvarGetLbGlobal(var));
2040 
2041  if ( (! negated && ! SCIPisEQ(scip, SCIPvarGetLbGlobal(var), 0.0)) || (negated && SCIPisEQ(scip, SCIPvarGetLbGlobal(var), 0.0)) )
2042  {
2043  /* the nonzeros are saved to later be inserted into the constant part (this is only done after all nonzeros of fixed variables have been
2044  * assembled, because we need to sort the constant nonzeros and loop over them, which we only want to do once and not once for each fixed
2045  * variable) */
2046  for (i = 0; i < consdata->nvarnonz[v]; i++)
2047  {
2048  savedcol[nfixednonz] = consdata->col[v][i];
2049  savedrow[nfixednonz] = consdata->row[v][i];
2050 
2051  /* this is the final value to add, we no longer have to remember from which variable this comes, minus because we have +A_i but -A_0 */
2052  if ( ! negated )
2053  savedval[nfixednonz] = consdata->val[v][i] * SCIPvarGetLbGlobal(var);
2054  else
2055  savedval[nfixednonz] = consdata->val[v][i]; /* if it is the negation of a variable fixed to zero, this variable is fixed to one */
2056 
2057  nfixednonz++;
2058  consdata->nnonz--;
2059  }
2060  }
2061  else
2062  {
2063  /* if the variable is fixed to zero, the nonzeros will just vanish, so we only reduce the number of nonzeros */
2064  consdata->nnonz -= consdata->nvarnonz[v];
2065  }
2066 
2067  /* free the memory of the corresponding entries in col/row/val */
2068  SCIPfreeBlockMemoryArrayNull(scip, &(consdata->val[v]), consdata->nvarnonz[v]);
2069  SCIPfreeBlockMemoryArrayNull(scip, &(consdata->row[v]), consdata->nvarnonz[v]);
2070  SCIPfreeBlockMemoryArrayNull(scip, &(consdata->col[v]), consdata->nvarnonz[v]);
2071 
2072  /* unlock variable */
2073  SCIP_CALL( unlockVar(scip, consdata, v) );
2074 
2075  /* as the variables don't need to be sorted, we just put the last variable into the empty spot and decrease sizes by one (at the end) */
2076  SCIP_CALL( SCIPreleaseVar(scip, &(consdata->vars[v])) );
2077  consdata->col[v] = consdata->col[consdata->nvars - 1];
2078  consdata->row[v] = consdata->row[consdata->nvars - 1];
2079  consdata->val[v] = consdata->val[consdata->nvars - 1];
2080  consdata->nvarnonz[v] = consdata->nvarnonz[consdata->nvars - 1];
2081  consdata->vars[v] = consdata->vars[consdata->nvars - 1];
2082  consdata->locks[v] = consdata->locks[consdata->nvars - 1];
2083  consdata->nvars--;
2084  v--; /* we need to check again if the variable we just shifted to this position also needs to be fixed */
2085  }
2086  else if ( aggregate && (SCIPvarGetStatus(var) == SCIP_VARSTATUS_AGGREGATED || SCIPvarGetStatus(var) == SCIP_VARSTATUS_MULTAGGR) )
2087  {
2088  SCIP_CALL( SCIPallocBufferArray(scip, &aggrvars, globalnvars) );
2089  SCIP_CALL( SCIPallocBufferArray(scip, &scalars, globalnvars) );
2090 
2091  /* this is how they should be initialized before calling SCIPgetProbvarLinearSum */
2092  if ( ! negated )
2093  {
2094  aggrvars[0] = consdata->vars[v];
2095  naggrvars = 1;
2096  constant = 0.0;
2097  scalars[0] = 1.0;
2098  }
2099  else
2100  {
2101  /* if this variable is the negation of var, than we look for a representation of 1.0 - var */
2102  aggrvars[0] = consdata->vars[v];
2103  naggrvars = 1;
2104  constant = 1.0;
2105  scalars[0] = -1.0;
2106  }
2107 
2108  /* get the variables this var was aggregated to */
2109  SCIP_CALL( SCIPgetProbvarLinearSum(scip, aggrvars, scalars, &naggrvars, globalnvars, &constant, &requiredsize, TRUE) );
2110  assert( requiredsize <= globalnvars ); /* requiredsize is the number of empty spots in aggrvars needed, globalnvars is the number
2111  * of spots we provided */
2112 
2113  /* Debugmessages for the (multi-)aggregation */
2114 #ifdef SCIP_DEBUG
2115  if ( SCIPvarGetStatus(consdata->vars[v]) == SCIP_VARSTATUS_AGGREGATED )
2116  SCIPdebugMsg(scip, "aggregating variable %s to ", SCIPvarGetName(var));
2117  else
2118  SCIPdebugMsg(scip, "multiaggregating variable %s to ", SCIPvarGetName(var));
2119  for (i = 0; i < naggrvars; i++)
2120  SCIPdebugMessagePrint(scip, "+ %g %s ", scalars[i], SCIPvarGetName(aggrvars[i]));
2121  SCIPdebugMessagePrint(scip, "+ %g.\n", constant);
2122 #endif
2123 
2124  /* add the nonzeros to the saved-arrays for the constant part, remove the nonzeros for the old variables and add them to the variables this variable
2125  * was (multi-)aggregated to */
2126  SCIP_CALL( multiaggrVar(scip, conss[c], v, aggrvars, scalars, naggrvars, constant, savedcol, savedrow, savedval, &nfixednonz, &vararraylength) );
2127  v--; /* we need to check again if the variable we just shifted to this position also needs to be fixed */
2128 
2129  SCIPfreeBufferArray(scip, &aggrvars);
2130  SCIPfreeBufferArray(scip, &scalars);
2131  }
2132  else if ( negated && (SCIPvarGetStatus(var) == SCIP_VARSTATUS_LOOSE || SCIPvarGetStatus(var) == SCIP_VARSTATUS_COLUMN) && aggregate )
2133  {
2134  /* if var1 is the negation of var2, then this is equivalent to it being aggregated to -var2 + 1 = 1 - var2 */
2135 
2136  SCIPdebugMsg(scip, "Changing variable %s to negation of variable <%s>!\n", SCIPvarGetName(consdata->vars[v]), SCIPvarGetName(var));
2137 
2138  scalar = -1.0;
2139 
2140  SCIP_CALL( multiaggrVar(scip, conss[c], v, &var, &scalar, 1, 1.0, savedcol, savedrow, savedval, &nfixednonz, &vararraylength) );
2141  v--; /* we need to check again if the variable we just shifted to this position also needs to be fixed */
2142  }
2143  }
2144 
2145  /* shrink the variable arrays if they were enlarged too much (or more vars were removed than added) */
2146  assert( consdata->nvars <= vararraylength );
2147  if ( consdata->nvars < vararraylength )
2148  {
2149  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->col, vararraylength, consdata->nvars) );
2150  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->row, vararraylength, consdata->nvars) );
2151  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->val, vararraylength, consdata->nvars) );
2152  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->nvarnonz, vararraylength, consdata->nvars) );
2153  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->vars, vararraylength, consdata->nvars) );
2154  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->locks, vararraylength, consdata->nvars) );
2155  }
2156 
2157  /* allocate the maximally needed memory for inserting the fixed variables into the constant part */
2158  arraylength = consdata->constnnonz + nfixednonz;
2159  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(consdata->constcol), consdata->constnnonz, arraylength) );
2160  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(consdata->constrow), consdata->constnnonz, arraylength) );
2161  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(consdata->constval), consdata->constnnonz, arraylength) );
2162 
2163  /* insert the fixed variables into the constant arrays, as we have +A_i but -A_0 we mutliply them by -1 */
2164  SCIP_CALL( SCIPsdpVarfixerMergeArrays(SCIPblkmem(scip), SCIPepsilon(scip), savedrow, savedcol, savedval, nfixednonz, FALSE, -1.0, consdata->constrow,
2165  consdata->constcol, consdata->constval, &(consdata->constnnonz), arraylength) );
2166 
2167  assert( consdata->constnnonz <= arraylength ); /* the allocated memory should always be sufficient */
2168 
2169  /* shrink the arrays if nonzeros could be combined */
2170  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(consdata->constcol), arraylength, consdata->constnnonz) );
2171  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(consdata->constrow), arraylength, consdata->constnnonz) );
2172  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(consdata->constval), arraylength, consdata->constnnonz) );
2173 
2174  /* free the saved arrays */
2175  SCIPfreeBufferArray(scip, &savedval);
2176  SCIPfreeBufferArray(scip, &savedrow);
2177  SCIPfreeBufferArray(scip, &savedcol);
2178 
2179  /* recompute sdpnnonz */
2180  consdata->nnonz = 0;
2181  for (v = 0; v < consdata->nvars; v++)
2182  consdata->nnonz += consdata->nvarnonz[v];
2183  }
2184 
2185  return SCIP_OKAY;
2186 }
2187 
2189 static
2190 SCIP_RETCODE enforceConstraint(
2191  SCIP* scip,
2192  SCIP_CONSHDLR* conshdlr,
2193  SCIP_CONS** conss,
2194  int nconss,
2195  SCIP_SOL* sol,
2196  SCIP_RESULT* result
2197  )
2198 {
2199  char cutname[SCIP_MAXSTRLEN];
2200  SCIP_CONSDATA* consdata;
2201  SCIP_CONSHDLRDATA* conshdlrdata;
2202  SCIP_Bool separated = FALSE;
2203  SCIP_ROW* row;
2204  SCIP_Bool infeasible;
2205  SCIP_Real lhs;
2206  SCIP_Real* coeff;
2207  SCIP_Real rhs;
2208  int nvars;
2209  int i;
2210  int j;
2211 #ifndef NDEBUG
2212  int snprintfreturn; /* used to check the return code of snprintf */
2213 #endif
2214 
2215  *result = SCIP_FEASIBLE;
2216 
2217  for (i = 0; i < nconss; ++i)
2218  {
2219  consdata = SCIPconsGetData(conss[i]);
2220  SCIP_CALL( SCIPconsSdpCheckSdpCons(scip, conss[i], sol, FALSE, result) );
2221  if ( *result == SCIP_FEASIBLE )
2222  continue;
2223 
2224  nvars = consdata->nvars;
2225 
2226  /* if the number of variables is 0, but the constraint is not feasible, we can cut off the node */
2227  if ( nvars == 0 )
2228  {
2229  *result = SCIP_CUTOFF;
2230  return SCIP_OKAY;
2231  }
2232 
2233  lhs = 0.0;
2234  coeff = NULL;
2235 
2236  SCIP_CALL( SCIPallocBufferArray(scip, &coeff, nvars) );
2237  SCIP_CALL( cutUsingEigenvector(scip, conss[i], sol, coeff, &lhs) );
2238 
2239  rhs = SCIPinfinity(scip);
2240  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2241 
2242 #ifndef NDEBUG
2243  snprintfreturn = SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "sepa_eig_sdp_%d", ++(conshdlrdata->neigveccuts));
2244  assert( snprintfreturn < SCIP_MAXSTRLEN ); /* check whether the name fits into the string */
2245 #else
2246  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "sepa_eig_sdp_%d", ++(conshdlrdata->neigveccuts));
2247 #endif
2248 
2249 #if ( SCIP_VERSION >= 700 || (SCIP_VERSION >= 602 && SCIP_SUBVERSION > 0) )
2250  SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, cutname , lhs, rhs, FALSE, FALSE, TRUE) );
2251 #else
2252  SCIP_CALL( SCIPcreateEmptyRowCons(scip, &row, conshdlr, cutname , lhs, rhs, FALSE, FALSE, TRUE) );
2253 #endif
2254  SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
2255 
2256  for (j = 0; j < nvars; ++j)
2257  {
2258  SCIP_CALL( SCIPaddVarToRow(scip, row, consdata->vars[j], coeff[j]) );
2259  }
2260 
2261  SCIP_CALL( SCIPflushRowExtensions(scip, row) );
2262 
2263 #ifdef SCIP_MORE_DEBUG
2264  SCIPinfoMessage(scip, NULL, "Added cut %s: ", cutname);
2265  SCIPinfoMessage(scip, NULL, "%f <= ", lhs);
2266  for (j = 0; j < nvars; j++)
2267  SCIPinfoMessage(scip, NULL, "+ (%f)*%s", coeff[j], SCIPvarGetName(consdata->vars[j]));
2268  SCIPinfoMessage(scip, NULL, "\n");
2269 #endif
2270 
2271  SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
2272 
2273  if ( infeasible )
2274  {
2275  *result = SCIP_CUTOFF;
2276 
2277  SCIP_CALL( SCIPreleaseRow(scip, &row) );
2278  SCIPfreeBufferArray(scip, &coeff);
2279 
2280  return SCIP_OKAY;
2281  }
2282  else
2283  {
2284  SCIP_CALL( SCIPaddPoolCut(scip, row) );
2285 
2286  SCIP_CALL( SCIPresetConsAge(scip, conss[i]) );
2287  *result = SCIP_SEPARATED;
2288  separated = TRUE;
2289  }
2290  SCIP_CALL( SCIPreleaseRow(scip, &row) );
2291  SCIPfreeBufferArray(scip, &coeff);
2292  }
2293 
2294  if ( separated )
2295  *result = SCIP_SEPARATED;
2296 
2297  return SCIP_OKAY;
2298 }
2299 
2300 
2302 static
2303 SCIP_DECL_QUADCONSUPGD(consQuadConsUpgdSdp)
2304 {
2305  char name[SCIP_MAXSTRLEN];
2306  SCIP_CONSHDLR* conshdlr;
2307  SCIP_CONSHDLRDATA* conshdlrdata;
2308  SCIP_CONS* lincons;
2309  SCIP_VAR** linconsvars;
2310  SCIP_Real* linconsvals;
2311  SCIP_VAR** linvarsterms;
2312  SCIP_Real* linvalsterms;
2313  SCIP_QUADVARTERM* quadvarterms;
2314  SCIP_BILINTERM* bilinterms;
2315  int nlinvarterms;
2316  int nquadvarterms;
2317  int nbilinterms;
2318  int nlinconsterms;
2319  int j;
2320 
2321  assert( scip != NULL );
2322  assert( cons != NULL );
2323  assert( nupgdconss != NULL );
2324  assert( upgdconss != NULL );
2325 
2326  *nupgdconss = 0;
2327 
2328  /* do not upgrade modifiable/sticking at node constraints */
2329  if ( SCIPconsIsModifiable(cons) || SCIPconsIsStickingAtNode(cons) )
2330  return SCIP_OKAY;
2331 
2332  /* do not run in sub-SCIPs to avoid recursive reformulations due to rank 1 constraints */
2333  if ( SCIPgetSubscipDepth(scip) > 0 )
2334  return SCIP_OKAY;
2335 
2336  /* make sure there is enough space to store the replacing constraints */
2337  if ( upgdconsssize < 1 )
2338  {
2339  *nupgdconss = -1;
2340  return SCIP_OKAY;
2341  }
2342 
2343  conshdlr = SCIPfindConshdlr(scip, CONSHDLRRANK1_NAME);
2344  if ( conshdlr == NULL )
2345  {
2346  SCIPerrorMessage("rank 1 SDP constraint handler not found\n");
2347  return SCIP_PLUGINNOTFOUND;
2348  }
2349  assert( conshdlr != NULL );
2350 
2351  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2352  assert( conshdlrdata != NULL );
2353 
2354  /* check whether upgrading should be performed */
2355  if ( ! conshdlrdata->sdpconshdlrdata->upgradquadconss )
2356  return SCIP_OKAY;
2357 
2358  /* we have to collect all variables appearing in quadratic constraints first */
2359  if ( conshdlrdata->sdpconshdlrdata->quadconsvars == NULL )
2360  {
2361  SCIP_CONSHDLR* quadconshdlr;
2362  SCIP_CONS** conss;
2363  int nconss;
2364  int nvars;
2365  int c;
2366  int i;
2367  int nsdpvars = 0;
2368 
2369  int** cols;
2370  int** rows;
2371  SCIP_Real** vals;
2372  SCIP_VAR** vars;
2373  int* nvarnonz;
2374  int nnonz;
2375  int nvarscnt;
2376  int constcol = 0;
2377  int constrow = 0;
2378  SCIP_Real constval = -1.0;
2379 
2380  /* todo: The arrays quadconsidx and quadconsvars are needed to check if variables have already been seen in a
2381  quadratic constraint. This could be replaced with a hashmap. */
2382  nvars = SCIPgetNTotalVars(scip);
2383  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &conshdlrdata->sdpconshdlrdata->quadconsidx, nvars) );
2384  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &conshdlrdata->sdpconshdlrdata->quadconsvars, nvars) );
2385  conshdlrdata->sdpconshdlrdata->nquadconsidx = nvars;
2386  for (j = 0; j < nvars; ++j)
2387  conshdlrdata->sdpconshdlrdata->quadconsidx[j] = -1;
2388 
2389  quadconshdlr = SCIPfindConshdlr(scip, "quadratic");
2390  if ( quadconshdlr == NULL )
2391  {
2392  SCIPerrorMessage("Quadratic constraint handler not found\n");
2393  return SCIP_PLUGINNOTFOUND;
2394  }
2395  assert( quadconshdlr != NULL );
2396 
2397  conss = SCIPconshdlrGetConss(quadconshdlr);
2398  nconss = SCIPconshdlrGetNConss(quadconshdlr);
2399 
2400  /* Do not perform upgrade, if there are too many quadratic constraints present. */
2401  if ( nconss > conshdlrdata->sdpconshdlrdata->maxnvarsquadupgd )
2402  {
2403  SCIPdebugMsg(scip, "There are %d many quadratic constraints present in the problem, thus do not upgrade quadratic constraints to an SDPrank1 constraint\n", nconss);
2404  SCIPfreeBlockMemoryArray(scip, &conshdlrdata->sdpconshdlrdata->quadconsvars, nvars);
2405  SCIPfreeBlockMemoryArray(scip, &conshdlrdata->sdpconshdlrdata->quadconsidx, nvars);
2406  return SCIP_OKAY;
2407  }
2408 
2409  for (c = 0; c < nconss; ++c)
2410  {
2411  assert( conss[c] != NULL );
2412 #ifdef SCIP_MORE_DEBUG
2413  SCIPinfoMessage(scip, NULL, "Found quadratic constraint to upgrade:\n");
2414  SCIP_CALL( SCIPprintCons(scip, conss[c], NULL) );
2415  SCIPinfoMessage(scip, NULL, "\n");
2416 #endif
2417  nquadvarterms = SCIPgetNQuadVarTermsQuadratic(scip, conss[c]);
2418  quadvarterms = SCIPgetQuadVarTermsQuadratic(scip, conss[c]);
2419 
2420  for (i = 0; i < nquadvarterms; ++i)
2421  {
2422  SCIP_VAR* var;
2423  int idx;
2424 
2425  assert( quadvarterms != NULL );
2426  var = quadvarterms[i].var;
2427  idx = SCIPvarGetIndex(var);
2428  assert( 0 <= idx && idx < nvars );
2429  if ( conshdlrdata->sdpconshdlrdata->quadconsidx[idx] < 0 )
2430  {
2431  conshdlrdata->sdpconshdlrdata->quadconsvars[nsdpvars] = var;
2432  conshdlrdata->sdpconshdlrdata->quadconsidx[idx] = nsdpvars++;
2433  }
2434  }
2435 
2436  nbilinterms = SCIPgetNBilinTermsQuadratic(scip, conss[c]);
2437  bilinterms = SCIPgetBilinTermsQuadratic(scip, conss[c]);
2438 
2439  for (i = 0; i < nbilinterms; ++i)
2440  {
2441  SCIP_VAR* var;
2442  int idx;
2443 
2444  assert( bilinterms != NULL );
2445  var = bilinterms[i].var1;
2446  idx = SCIPvarGetIndex(var);
2447  assert( 0 <= idx && idx < nvars );
2448  if ( conshdlrdata->sdpconshdlrdata->quadconsidx[idx] < 0 )
2449  {
2450  conshdlrdata->sdpconshdlrdata->quadconsvars[nsdpvars] = var;
2451  conshdlrdata->sdpconshdlrdata->quadconsidx[idx] = nsdpvars++;
2452  }
2453 
2454  var = bilinterms[i].var2;
2455  idx = SCIPvarGetIndex(var);
2456  assert( 0 <= idx && idx < nvars );
2457  if ( conshdlrdata->sdpconshdlrdata->quadconsidx[idx] < 0 )
2458  {
2459  conshdlrdata->sdpconshdlrdata->quadconsvars[nsdpvars] = var;
2460  conshdlrdata->sdpconshdlrdata->quadconsidx[idx] = nsdpvars++;
2461  }
2462  }
2463  }
2464 
2465  /* do not perform upgrade, if there are too many variables in the quadratic constraints, since we need sdpvars *
2466  sdpvars many variables for the (dual) SDPrank1 constraint */
2467  if ( nsdpvars > conshdlrdata->sdpconshdlrdata->maxnvarsquadupgd )
2468  {
2469  SCIPdebugMsg(scip, "There are %d many variables present in the quadratic constraints, thus do not upgrade quadratic constraints to an SDPrank1 constraint\n", nsdpvars);
2470  SCIPfreeBlockMemoryArray(scip, &conshdlrdata->sdpconshdlrdata->quadconsvars, nvars);
2471  SCIPfreeBlockMemoryArray(scip, &conshdlrdata->sdpconshdlrdata->quadconsidx, nvars);
2472  return SCIP_OKAY;
2473  }
2474 
2475  /* create bilinear variables */
2476  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &conshdlrdata->sdpconshdlrdata->X, nsdpvars) );
2477  conshdlrdata->sdpconshdlrdata->nsdpvars = nsdpvars;
2478 
2479  for (i = 0; i < nsdpvars; ++i)
2480  {
2481  SCIP_Real lb1;
2482  SCIP_Real ub1;
2483  SCIP_VAR* var1;
2484 
2485  var1 = conshdlrdata->sdpconshdlrdata->quadconsvars[i];
2486  assert( var1 != NULL );
2487  lb1 = SCIPvarGetLbGlobal(var1);
2488  ub1 = SCIPvarGetUbGlobal(var1);
2489 
2490  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &conshdlrdata->sdpconshdlrdata->X[i], nsdpvars) );
2491 
2492  for (j = 0; j <= i; ++j)
2493  {
2494  SCIP_VARTYPE vartype;
2495  SCIP_VAR* var2;
2496  SCIP_Real lb2;
2497  SCIP_Real ub2;
2498  SCIP_Real lb;
2499  SCIP_Real ub;
2500 
2501  var2 = conshdlrdata->sdpconshdlrdata->quadconsvars[j];
2502  assert( var2 != NULL );
2503  lb2 = SCIPvarGetLbGlobal(var2);
2504  ub2 = SCIPvarGetUbGlobal(var2);
2505 
2506  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "X%d#%d", i, j);
2507 
2508  lb = MIN3(lb1 * lb2, lb1 * ub2, ub1 * lb2);
2509  lb = MIN(lb, ub1 * ub2);
2510  ub = MAX3(lb1 * lb2, lb1 * ub2, ub1 * lb2);
2511  ub = MAX(ub, ub1 * ub2);
2512 
2513  if ( SCIPvarIsBinary(var1) && SCIPvarIsBinary(var2) )
2514  vartype = SCIP_VARTYPE_BINARY;
2515  else if ( SCIPvarIsIntegral(var1) && SCIPvarIsIntegral(var2) )
2516  vartype = SCIP_VARTYPE_INTEGER;
2517  else
2518  vartype = SCIP_VARTYPE_CONTINUOUS;
2519 
2520  SCIP_CALL( SCIPcreateVarBasic(scip, &(conshdlrdata->sdpconshdlrdata->X[i][j]), name, lb, ub, 0.0, vartype) );
2521  SCIP_CALL( SCIPaddVar(scip, conshdlrdata->sdpconshdlrdata->X[i][j]) );
2522  }
2523  }
2524 
2525  /* fill SDP data */
2526  nnonz = nsdpvars + nsdpvars * (nsdpvars + 1) / 2;
2527  SCIP_CALL( SCIPallocBufferArray(scip, &cols, nnonz) );
2528  SCIP_CALL( SCIPallocBufferArray(scip, &rows, nnonz) );
2529  SCIP_CALL( SCIPallocBufferArray(scip, &vals, nnonz) );
2530  SCIP_CALL( SCIPallocBufferArray(scip, &vars, nnonz) );
2531  SCIP_CALL( SCIPallocBufferArray(scip, &nvarnonz, nnonz) );
2532 
2533  /* first the terms for the original variables */
2534  for (j = 0; j < nsdpvars; ++j)
2535  {
2536  SCIP_CALL( SCIPallocBufferArray(scip, &cols[j], 1) );
2537  SCIP_CALL( SCIPallocBufferArray(scip, &rows[j], 1) );
2538  SCIP_CALL( SCIPallocBufferArray(scip, &vals[j], 1) );
2539  nvarnonz[j] = 1;
2540  cols[j][0] = 0;
2541  rows[j][0] = 1 + j;
2542  vals[j][0] = 1.0;
2543  vars[j] = conshdlrdata->sdpconshdlrdata->quadconsvars[j];
2544  }
2545 
2546  /* now the terms for the bilinear terms */
2547  nvarscnt = nsdpvars;
2548  for (i = 0; i < nsdpvars; ++i)
2549  {
2550  for (j = 0; j <= i; ++j)
2551  {
2552  SCIP_CALL( SCIPallocBufferArray(scip, &cols[nvarscnt], 1) );
2553  SCIP_CALL( SCIPallocBufferArray(scip, &rows[nvarscnt], 1) );
2554  SCIP_CALL( SCIPallocBufferArray(scip, &vals[nvarscnt], 1) );
2555  nvarnonz[nvarscnt] = 1;
2556  cols[nvarscnt][0] = 1 + j;
2557  rows[nvarscnt][0] = 1 + i;
2558  vals[nvarscnt][0] = 1.0;
2559  vars[nvarscnt] = conshdlrdata->sdpconshdlrdata->X[i][j];
2560  ++nvarscnt;
2561  }
2562  }
2563  assert( nvarscnt == nsdpvars + nsdpvars * (nsdpvars + 1)/2 );
2564 
2565  /* create corresponding rank 1 SDP constraint */
2566  if ( conshdlrdata->sdpconshdlrdata->upgradekeepquad )
2567  {
2568  SCIP_CALL( SCIPcreateConsSdp(scip, &conshdlrdata->sdpconshdlrdata->sdpcons, "QuadraticSDPcons", nvarscnt, nvarscnt, 1 + nsdpvars, nvarnonz,
2569  cols, rows, vals, vars, 1, &constcol, &constrow, &constval) );
2570  SCIP_CALL( SCIPaddCons(scip, conshdlrdata->sdpconshdlrdata->sdpcons) );
2571  }
2572  else
2573  {
2574  SCIP_CALL( SCIPcreateConsSdpRank1(scip, &conshdlrdata->sdpconshdlrdata->sdpcons, "QuadraticSDPrank1cons", nvarscnt, nvarscnt, 1 + nsdpvars, nvarnonz,
2575  cols, rows, vals, vars, 1, &constcol, &constrow, &constval) );
2576  SCIP_CALL( SCIPaddCons(scip, conshdlrdata->sdpconshdlrdata->sdpcons) );
2577  }
2578 
2579 #ifdef SCIP_MORE_DEBUG
2580  SCIPinfoMessage(scip, NULL, "In upgrade of quadratic constraint the following SDPrank1 constraint has been added:\n");
2581  SCIP_CALL( SCIPprintCons(scip, conshdlrdata->sdpconshdlrdata->sdpcons, NULL) );
2582  SCIPinfoMessage(scip, NULL, "\n");
2583 #endif
2584 
2585  /* free local memory */
2586  for (j = nvarscnt - 1; j >= 0; --j)
2587  {
2588  SCIPfreeBufferArray(scip, &vals[j]);
2589  SCIPfreeBufferArray(scip, &rows[j]);
2590  SCIPfreeBufferArray(scip, &cols[j]);
2591  }
2592  SCIPfreeBufferArray(scip, &nvarnonz);
2593  SCIPfreeBufferArray(scip, &vars);
2594  SCIPfreeBufferArray(scip, &vals);
2595  SCIPfreeBufferArray(scip, &rows);
2596  SCIPfreeBufferArray(scip, &cols);
2597  }
2598 
2599  if ( ! conshdlrdata->sdpconshdlrdata->upgradekeepquad )
2600  {
2601  int cnt = 0;
2602 
2603  /* create linear constraint for quadratic constraint */
2604  nlinvarterms = SCIPgetNLinearVarsQuadratic(scip, cons);
2605  linvarsterms = SCIPgetLinearVarsQuadratic(scip, cons);
2606  linvalsterms = SCIPgetCoefsLinearVarsQuadratic(scip, cons);
2607  nquadvarterms = SCIPgetNQuadVarTermsQuadratic(scip, cons);
2608  quadvarterms = SCIPgetQuadVarTermsQuadratic(scip, cons);
2609  nbilinterms = SCIPgetNBilinTermsQuadratic(scip, cons);
2610  bilinterms = SCIPgetBilinTermsQuadratic(scip, cons);
2611 
2612  /* a quadvarterm consists of a variable x and two coefficients, one for the linear term x and one for the quadratic
2613  term x^2, where at least one of the two coefficients is nonzero */
2614  nlinconsterms = nlinvarterms + 2 * nquadvarterms + nbilinterms;
2615  SCIP_CALL( SCIPallocBufferArray(scip, &linconsvars, nlinconsterms) );
2616  SCIP_CALL( SCIPallocBufferArray(scip, &linconsvals, nlinconsterms) );
2617 
2618  /* fill in constraint */
2619  for (j = 0; j < nlinvarterms; ++j)
2620  {
2621  linconsvals[cnt] = linvalsterms[j];
2622  linconsvars[cnt] = linvarsterms[j];
2623  assert( linconsvars[cnt] != NULL );
2624  ++cnt;
2625  }
2626  assert( cnt == nlinvarterms );
2627  for (j = 0; j < nquadvarterms; ++j)
2628  {
2629  int idx;
2630 
2631  idx = SCIPvarGetIndex(quadvarterms[j].var);
2632  idx = conshdlrdata->sdpconshdlrdata->quadconsidx[idx];
2633  assert( 0 <= idx && idx < conshdlrdata->sdpconshdlrdata->nsdpvars );
2634 
2635  /* add coefficient for linear term corresponding to the current variable (may be zero) */
2636  if ( ! SCIPisZero(scip, quadvarterms[j].lincoef) )
2637  {
2638  linconsvals[cnt] = quadvarterms[j].lincoef;
2639  linconsvars[cnt] = quadvarterms[j].var;
2640  assert( linconsvars[cnt] != NULL );
2641  ++cnt;
2642  }
2643 
2644  /* add coefficient for quadratic term corresponding to the current variable (may be zero) */
2645  if ( ! SCIPisZero(scip, quadvarterms[j].sqrcoef) )
2646  {
2647  linconsvals[cnt] = quadvarterms[j].sqrcoef;
2648  linconsvars[cnt] = conshdlrdata->sdpconshdlrdata->X[idx][idx];
2649  assert( linconsvars[cnt] != NULL );
2650  ++cnt;
2651  }
2652 
2653  SCIPdebugMsg(scip, "New variable %s corresponds to squared original variable %s\n", SCIPvarGetName(conshdlrdata->sdpconshdlrdata->X[idx][idx]), SCIPvarGetName(quadvarterms[j].var));
2654  }
2655  assert( cnt <= nlinvarterms + 2 * nquadvarterms );
2656 
2657  for (j = 0; j < nbilinterms; ++j)
2658  {
2659  int idx1;
2660  int idx2;
2661 
2662  idx1 = SCIPvarGetIndex(bilinterms[j].var1);
2663  idx1 = conshdlrdata->sdpconshdlrdata->quadconsidx[idx1];
2664  assert( 0 <= idx1 && idx1 < conshdlrdata->sdpconshdlrdata->nsdpvars );
2665 
2666  idx2 = SCIPvarGetIndex(bilinterms[j].var2);
2667  idx2 = conshdlrdata->sdpconshdlrdata->quadconsidx[idx2];
2668  assert( 0 <= idx2 && idx2 < conshdlrdata->sdpconshdlrdata->nsdpvars );
2669 
2670  if ( idx2 > idx1 )
2671  SCIPswapInts(&idx1, &idx2);
2672 
2673  linconsvals[cnt] = bilinterms[j].coef;
2674  linconsvars[cnt] = conshdlrdata->sdpconshdlrdata->X[idx1][idx2];
2675  assert( linconsvars[cnt] != NULL );
2676  ++cnt;
2677 
2678  SCIPdebugMsg(scip, "New variable %s corresponds to product of original variables %s and %s\n", SCIPvarGetName(conshdlrdata->sdpconshdlrdata->X[idx1][idx2]), SCIPvarGetName(bilinterms[j].var1), SCIPvarGetName(bilinterms[j].var2));
2679  }
2680  assert( cnt <= nlinvarterms + 2 * nquadvarterms + nbilinterms );
2681 
2682  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "lin_%s", SCIPconsGetName(cons));
2683  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, name, cnt, linconsvars, linconsvals, SCIPgetLhsQuadratic(scip, cons), SCIPgetRhsQuadratic(scip, cons),
2684  SCIPconsIsInitial(cons), SCIPconsIsSeparated(cons), SCIPconsIsEnforced(cons), SCIPconsIsChecked(cons), SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons),
2685  FALSE, SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons), FALSE) );
2686 
2687 #ifdef SCIP_MORE_DEBUG
2688  SCIPinfoMessage(scip, NULL, "In upgrade of quadratic constraint the following linear constraint has been added:\n");
2689  SCIP_CALL( SCIPprintCons(scip, lincons, NULL) );
2690  SCIPinfoMessage(scip, NULL, "\n");
2691 #endif
2692 
2693  /* fill in upgdconss - do not mention SDP constraint, since this has been added already */
2694  upgdconss[0] = lincons;
2695  *nupgdconss = 1;
2696 
2697  SCIPfreeBufferArray(scip, &linconsvals);
2698  SCIPfreeBufferArray(scip, &linconsvars);
2699  }
2700  else
2701  {
2702  /* todo: Check whether adding the linear constraints helps */
2703  *nupgdconss = 0; /* the original quadratic constraint should be kept in the problem */
2704  }
2705 
2706  return SCIP_OKAY;
2707 }
2708 
2709 
2710 /*
2711  * callbacks
2712  */
2713 
2715 static
2716 SCIP_DECL_CONSINITPRE(consInitpreSdp)
2717 {/*lint --e{715}*/
2718  SCIP_CONSHDLRDATA* conshdlrdata;
2719 
2720  assert( conshdlr != NULL );
2721 
2722  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2723  assert( conshdlrdata != NULL );
2724 
2725  conshdlrdata->neigveccuts = 0; /* this is used to give the eigenvector-cuts distinguishable names */
2726  conshdlrdata->ndiaggezerocuts = 0; /* this is used to give the diagGEzero-cuts distinguishable names */
2727  conshdlrdata->n1x1blocks = 0; /* this is used to give the lp constraints resulting from 1x1 sdp-blocks distinguishable names */
2728 
2729  return SCIP_OKAY;
2730 }
2731 
2733 static
2734 SCIP_DECL_CONSLOCK(consLockSdp)
2735 {/*lint --e{715}*/
2736  SCIP_Real* Aj;
2737  SCIP_CONSDATA* consdata;
2738  int nvars;
2739  int v;
2740 
2741  consdata = SCIPconsGetData(cons);
2742  assert( consdata != NULL );
2743  nvars = consdata->nvars;
2744 
2745  SCIPdebugMsg(scip, "locking method of <%s>.\n", SCIPconsGetName(cons));
2746 
2747  /* rank-1 constraints are always up- and down-locked */
2748  if ( consdata->rankone )
2749  {
2750  if ( consdata->locks == NULL )
2751  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->locks, nvars) );
2752 
2753  for (v = 0; v < consdata->nvars; ++v)
2754  {
2755  consdata->locks[v] = 0;
2756  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], locktype, nlockspos + nlocksneg, nlockspos + nlocksneg) );
2757  }
2758  return SCIP_OKAY;
2759  }
2760 
2761  /* if locks have not yet been computed */
2762  if ( consdata->locks == NULL )
2763  {
2764  SCIP_Real eigenvalue;
2765  int blocksize;
2766 
2767  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->locks, nvars) );
2768 
2769  blocksize = consdata->blocksize;
2770 
2771  SCIP_CALL( SCIPallocBufferArray(scip, &Aj, blocksize * blocksize) ); /*lint !e647*/
2772 
2773  for (v = 0; v < nvars; v++)
2774  {
2775  SCIP_CALL( SCIPconsSdpGetFullAj(scip, cons, v, Aj) );
2776  consdata->locks[v] = -2; /* unintitialized */
2777 
2778  /* compute the smallest eigenvalue */
2779  SCIP_CALL( SCIPlapackComputeIthEigenvalue(SCIPbuffer(scip), FALSE, blocksize, Aj, 1, &eigenvalue, NULL) );
2780  if ( SCIPisNegative(scip, eigenvalue) )
2781  {
2782  /* as the lowest eigenvalue is negative, the matrix is not positive semidefinite, so adding more of it can remove positive
2783  * semidefiniteness of the SDP-matrix */
2784  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], locktype, nlocksneg, nlockspos) );
2785  consdata->locks[v] = 1; /* up-lock */
2786  }
2787 
2788  /* if the smallest eigenvalue is already positive, we don't need to compute the biggest one */
2789  if ( SCIPisPositive(scip, eigenvalue) )
2790  {
2791  /* as an eigenvalue is positive, the matrix is not negative semidefinite, so substracting more of it can remove positive
2792  * semidefiniteness of the SDP-matrix */
2793  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], locktype, nlockspos, nlocksneg) );
2794  consdata->locks[v] = -1; /* down-lock */
2795  }
2796  else
2797  {
2798  /* compute the biggest eigenvalue */
2799  SCIP_CALL( SCIPlapackComputeIthEigenvalue(SCIPbuffer(scip), FALSE, blocksize, Aj, blocksize, &eigenvalue, NULL) );
2800  if ( SCIPisPositive(scip, eigenvalue) )
2801  {
2802  /* as the biggest eigenvalue is positive, the matrix is not negative semidefinite, so substracting more of it can remove positive
2803  * semidefiniteness of the SDP-matrix */
2804  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], locktype, nlockspos, nlocksneg) );
2805  if ( consdata->locks[v] == 1 )
2806  {
2807  consdata->locks[v] = 0; /* up- and down-lock */
2808  }
2809  else
2810  consdata->locks[v] = -1; /* down-lock */
2811  }
2812  }
2813  }
2814 
2815  SCIPfreeBufferArray(scip, &Aj);
2816  }
2817  else
2818  {
2819 #ifndef NDEBUG
2820  SCIP_CALL( checkVarsLocks(scip, cons) );
2821 #endif
2822  for (v = 0; v < nvars; v++)
2823  {
2824  if ( consdata->locks[v] == 1 ) /* up-lock */
2825  {
2826  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], locktype, nlocksneg, nlockspos) );
2827  }
2828  else if ( consdata->locks[v] == -1 ) /* down-lock */
2829  {
2830  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], locktype, nlockspos, nlocksneg) );
2831  }
2832  else if ( consdata->locks[v] == 0 ) /* up and down lock */
2833  {
2834  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[v], locktype, nlockspos + nlocksneg, nlockspos + nlocksneg) );
2835  }
2836  else
2837  assert( consdata->locks[v] == -2 );
2838  }
2839  }
2840 
2841  return SCIP_OKAY;
2842 }
2843 
2848 static
2849 SCIP_DECL_CONSEXIT(consExitSdp)
2850 {
2851  SCIP_CONSHDLRDATA* conshdlrdata;
2852 
2853  assert(scip != NULL);
2854 
2855  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2856  assert(conshdlrdata != NULL);
2857 
2858  /* reset parameter triedlinearconss */
2859  conshdlrdata->sdpconshdlrdata->triedlinearconss = FALSE;
2860 
2861  return SCIP_OKAY;
2862 }
2863 
2865 static
2866 SCIP_DECL_CONSEXITPRE(consExitpreSdp)
2867 {/*lint --e{715}*/
2868  SCIP_CONSHDLRDATA* conshdlrdata;
2869  int i;
2870  int j;
2871 
2872  assert( scip != NULL );
2873  assert( conshdlr != NULL );
2874 
2875  SCIPdebugMsg(scip, "Exitpre method of conshdlr <%s>.\n", SCIPconshdlrGetName(conshdlr));
2876 
2877  if ( conss == NULL )
2878  return SCIP_OKAY;
2879 
2880  SCIP_CALL( fixAndAggrVars(scip, conss, nconss, TRUE) );
2881 
2882  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2883  assert( conshdlrdata != NULL );
2884 
2885  SCIPfreeBlockMemoryArrayNull(scip, &conshdlrdata->sdpconshdlrdata->quadconsidx, conshdlrdata->sdpconshdlrdata->nquadconsidx);
2886  SCIPfreeBlockMemoryArrayNull(scip, &conshdlrdata->sdpconshdlrdata->quadconsvars, conshdlrdata->sdpconshdlrdata->nquadconsidx);
2887  if ( conshdlrdata->sdpconshdlrdata->X != NULL )
2888  {
2889  SCIPdebugMsg(scip, "Releasing additional variables from upgrading method\n");
2890  for (i = 0; i < conshdlrdata->sdpconshdlrdata->nsdpvars; ++i)
2891  {
2892  for (j = 0; j <= i; ++j)
2893  {
2894  SCIP_CALL( SCIPreleaseVar(scip, &(conshdlrdata->sdpconshdlrdata->X[i][j])) );
2895  }
2896  SCIPfreeBlockMemoryArray(scip, &conshdlrdata->sdpconshdlrdata->X[i], conshdlrdata->sdpconshdlrdata->nsdpvars);
2897  }
2898  SCIPfreeBlockMemoryArray(scip, &conshdlrdata->sdpconshdlrdata->X, conshdlrdata->sdpconshdlrdata->nsdpvars);
2899  }
2900  if ( conshdlrdata->sdpconshdlrdata->sdpcons != NULL )
2901  {
2902  SCIPdebugMsg(scip, "Releasing constraint %s from upgrading method\n", SCIPconsGetName(conshdlrdata->sdpconshdlrdata->sdpcons) );
2903  SCIP_CALL( SCIPreleaseCons(scip, &conshdlrdata->sdpconshdlrdata->sdpcons) );
2904  }
2905 
2906  /* TODO: test if branching and/or separation of Chen et al. can be applied */
2907 
2908  return SCIP_OKAY;
2909 }
2910 
2915 static
2916 SCIP_DECL_CONSINITSOL(consInitsolSdp)
2917 {/*lint --e{715}*/
2918  SCIP_CONSHDLRDATA* conshdlrdata;
2919  int c;
2920  int i;
2921 
2922  assert( scip != NULL );
2923  assert( conshdlr != NULL );
2924 
2925  if ( conss == NULL )
2926  return SCIP_OKAY;
2927 
2928  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2929  assert( conshdlrdata != NULL );
2930 
2931  if ( SCIPgetSubscipDepth(scip) > 0 || ! conshdlrdata->sdpconshdlrdata->quadconsrank1 )
2932  return SCIP_OKAY;
2933 
2934  for (c = 0; c < nconss; ++c)
2935  {
2936  SCIP_CONSDATA* consdata;
2937 
2938  consdata = SCIPconsGetData(conss[c]);
2939  assert( consdata != NULL );
2940 
2941  consdata->maxevsubmat[0] = -1;
2942  consdata->maxevsubmat[1] = -1;
2943 
2944  /* For each constraint, if it should be rank one, add all quadratic constraints given by the 2x2 principal
2945  * minors. */
2946  if ( consdata->rankone && ! consdata->addedquadcons )
2947  {
2948  SCIP_VAR** quadvars1;
2949  SCIP_VAR** quadvars2;
2950  SCIP_VAR** linvars;
2951  SCIP_CONS* quadcons;
2952  SCIP_Real* lincoefs;
2953  SCIP_Real* quadcoefs;
2954  SCIP_Real* constmatrix;
2955  SCIP_Real** matrixAk;
2956  SCIP_Real lhs;
2957  SCIP_Real aiik;
2958  SCIP_Real ajjk;
2959  SCIP_Real aijk;
2960  SCIP_Real aiil;
2961  SCIP_Real ajjl;
2962  SCIP_Real aijl;
2963  SCIP_Real cii;
2964  SCIP_Real cjj;
2965  SCIP_Real cij;
2966  char name[SCIP_MAXSTRLEN];
2967  int* nnonzvars;
2968  int** nonzvars;
2969  int j;
2970  int k;
2971  int l;
2972  int blocksize;
2973  int varind1;
2974  int varind2;
2975 
2976  blocksize = consdata->blocksize;
2977 
2978  SCIP_CALL( SCIPallocBufferArray(scip, &constmatrix, (blocksize * (blocksize + 1)) / 2) ); /*lint !e647*/
2979  SCIP_CALL( SCIPconsSdpGetLowerTriangConstMatrix(scip, conss[c], constmatrix) );
2980 
2981  SCIP_CALL( SCIPallocBufferArray(scip, &quadvars1, consdata->nvars * consdata->nvars) );
2982  SCIP_CALL( SCIPallocBufferArray(scip, &quadvars2, consdata->nvars * consdata->nvars) );
2983  SCIP_CALL( SCIPallocBufferArray(scip, &linvars, consdata->nvars) );
2984  SCIP_CALL( SCIPallocBufferArray(scip, &quadcoefs, consdata->nvars * consdata->nvars) );
2985  SCIP_CALL( SCIPallocBufferArray(scip, &lincoefs, consdata->nvars) );
2986  SCIP_CALL( SCIPallocBufferArray(scip, &matrixAk, consdata->nvars) );
2987 
2988  for (i = 0; i < consdata->nvars; ++i)
2989  {
2990  SCIP_CALL( SCIPallocBufferArray(scip, &matrixAk[i], blocksize * blocksize) );
2991  SCIP_CALL( SCIPconsSdpGetFullAj(scip, conss[c], i, matrixAk[i]) );
2992  }
2993 
2994  SCIP_CALL( SCIPallocBufferArray(scip, &nnonzvars, (blocksize * (blocksize + 1)) / 2) );
2995  SCIP_CALL( SCIPallocBufferArray(scip, &nonzvars, (blocksize * (blocksize + 1)) / 2) );
2996 
2997  for (i = 0; i < blocksize; ++i)
2998  {
2999  for (j = 0; j <= i; ++j)
3000  {
3001  int varcnt = 0;
3002 
3003  SCIP_CALL( SCIPallocBufferArray(scip, &nonzvars[SCIPconsSdpCompLowerTriangPos(i,j)], consdata->nvars) );
3004 
3005  for (k = 0; k < consdata->nvars; ++k)
3006  {
3007  if ( ! SCIPisZero(scip, matrixAk[k][i * blocksize + j]) || ! SCIPisZero(scip, matrixAk[k][i * blocksize + i]) || ! SCIPisZero(scip, matrixAk[k][j * blocksize + j]) )
3008  {
3009  nonzvars[SCIPconsSdpCompLowerTriangPos(i,j)][varcnt] = k;
3010  varcnt++;
3011  }
3012  }
3013  nnonzvars[SCIPconsSdpCompLowerTriangPos(i,j)] = varcnt;
3014  }
3015  }
3016 
3017  for (i = 0; i < blocksize; ++i)
3018  {
3019  for (j = 0; j < i; ++j)
3020  {
3021  int lincnt = 0;
3022  int quadcnt = 0;
3023 
3024  cii = constmatrix[SCIPconsSdpCompLowerTriangPos(i,i)];
3025  cjj = constmatrix[SCIPconsSdpCompLowerTriangPos(j,j)];
3026  cij = constmatrix[SCIPconsSdpCompLowerTriangPos(i,j)];
3027 
3028  for (k = 0; k < nnonzvars[SCIPconsSdpCompLowerTriangPos(i,j)]; ++k)
3029  {
3030  varind1 = nonzvars[SCIPconsSdpCompLowerTriangPos(i,j)][k];
3031  ajjk = matrixAk[varind1][j * consdata->blocksize + j];
3032  aiik = matrixAk[varind1][i * consdata->blocksize + i];
3033  aijk = matrixAk[varind1][j * consdata->blocksize + i];
3034 
3035  if ( ! SCIPisZero(scip, -cii * ajjk - cjj * aiik + cij * aijk) )
3036  {
3037  linvars[lincnt] = consdata->vars[varind1];
3038  lincoefs[lincnt] = -cii * ajjk - cjj * aiik + cij * aijk;
3039  ++lincnt;
3040  }
3041 
3042  for (l = 0; l < k; ++l)
3043  {
3044  varind2 = nonzvars[SCIPconsSdpCompLowerTriangPos(i,j)][l];
3045  ajjl = matrixAk[varind2][j * consdata->blocksize + j];
3046  aiil = matrixAk[varind2][i * consdata->blocksize + i];
3047  aijl = matrixAk[varind2][j * consdata->blocksize + i];
3048 
3049  if ( ! SCIPisZero(scip, aiik * ajjl + ajjk * aiil - 2 * aijk * aijl) )
3050  {
3051  quadvars1[quadcnt] = consdata->vars[varind1];
3052  quadvars2[quadcnt] = consdata->vars[varind2];
3053  quadcoefs[quadcnt] = aiik * ajjl + ajjk * aiil - 2 * aijk * aijl;
3054  ++quadcnt;
3055  }
3056  }
3057 
3058  /* case l == k needs special treatment */
3059  if ( ! SCIPisZero(scip, aiik * ajjk - aijk * aijk) )
3060  {
3061  quadvars1[quadcnt] = consdata->vars[varind1];
3062  quadvars2[quadcnt] = consdata->vars[varind1];
3063  quadcoefs[quadcnt] = aiik * ajjk - aijk * aijk;
3064  ++quadcnt;
3065  }
3066  }
3067  assert( quadcnt <= consdata->nvars * consdata->nvars );
3068  assert( lincnt <= consdata->nvars );
3069 
3070  lhs = cij * cij - cii * cjj;
3071 
3072  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "quadcons#%d#%d#%d", i, j, c);
3073 
3074  /* create quadratic constraint */
3075  SCIP_CALL( SCIPcreateConsQuadratic(scip, &quadcons, name, lincnt, linvars, lincoefs, quadcnt, quadvars1, quadvars2, quadcoefs, lhs, lhs,
3076  FALSE, /* initial */
3077  TRUE, /* separate */
3078  TRUE, /* enforce */
3079  TRUE, /* check */
3080  TRUE, /* propagate */
3081  FALSE, /* local */
3082  FALSE, /* modifiable */
3083  FALSE, /* dynamic */
3084  TRUE) ); /* removable */
3085 
3086 #ifdef SCIP_MORE_DEBUG
3087  SCIP_CALL( SCIPprintCons(scip, quadcons, NULL) );
3088  SCIPinfoMessage(scip, NULL, "\n");
3089 #endif
3090 
3091  SCIP_CALL( SCIPaddCons(scip, quadcons) );
3092  SCIP_CALL( SCIPreleaseCons(scip, &quadcons) );
3093  }
3094  }
3095 
3096  for (i = 0; i < blocksize; ++i)
3097  {
3098  for (j = 0; j <= i; ++j)
3099  SCIPfreeBufferArray(scip, &nonzvars[SCIPconsSdpCompLowerTriangPos(i,j)]);
3100  }
3101 
3102  SCIPfreeBufferArray(scip, &nonzvars);
3103  SCIPfreeBufferArray(scip, &nnonzvars);
3104 
3105  for (i = 0; i < consdata->nvars; ++i)
3106  SCIPfreeBufferArray(scip, &matrixAk[i]);
3107 
3108  SCIPfreeBufferArray(scip, &matrixAk);
3109  SCIPfreeBufferArray(scip, &lincoefs);
3110  SCIPfreeBufferArray(scip, &quadcoefs);
3111  SCIPfreeBufferArray(scip, &linvars);
3112  SCIPfreeBufferArray(scip, &quadvars2);
3113  SCIPfreeBufferArray(scip, &quadvars1);
3114  SCIPfreeBufferArray(scip, &constmatrix);
3115  }
3116  consdata->addedquadcons = TRUE;
3117  }
3118  return SCIP_OKAY;
3119 }
3120 
3122 static
3123 SCIP_DECL_CONSPRESOL(consPresolSdp)
3124 {/*lint --e{715}*/
3125  SCIP_CONSHDLRDATA* conshdlrdata;
3126 
3127  assert( conshdlr != NULL );
3128  assert( result != NULL );
3129 
3130  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3131  assert( conshdlrdata != NULL );
3132 
3133  *result = SCIP_DIDNOTRUN;
3134 
3135  if ( nrounds == 0 )
3136  {
3137  int noldaddconss;
3138  int nolddelconss;
3139  int noldchgbds;
3140 
3141  *result = SCIP_DIDNOTFIND;
3142 
3143  noldaddconss = *naddconss;
3144  nolddelconss = *ndelconss;
3145  noldchgbds = *nchgbds;
3146  SCIP_CALL( move_1x1_blocks_to_lp(scip, conshdlr, conss, nconss, naddconss, ndelconss, nchgbds, result) );
3147  if ( *result != SCIP_CUTOFF && (noldaddconss != *naddconss || nolddelconss != *ndelconss || noldchgbds != *nchgbds) )
3148  *result = SCIP_SUCCESS;
3149 
3150  /* In the following, we add linear constraints to be propagated. This is needed only once. We assume that this is
3151  * only necessary in the main SCIP instance. */
3152  if ( SCIPgetSubscipDepth(scip) == 0 && ! conshdlrdata->sdpconshdlrdata->triedlinearconss )
3153  {
3154  conshdlrdata->sdpconshdlrdata->triedlinearconss = TRUE;
3155  if ( *result != SCIP_CUTOFF && conshdlrdata->sdpconshdlrdata->diaggezerocuts )
3156  {
3157  noldaddconss = *naddconss;
3158  noldchgbds = *nchgbds;
3159  SCIP_CALL( diagGEzero(scip, conshdlr, conss, nconss, naddconss, nchgbds, result) );
3160  SCIPdebugMsg(scip, "Diagonal entries: added %d cuts and changed %d bounds.\n", *naddconss - noldaddconss, *nchgbds - noldchgbds);
3161  if ( *result != SCIP_CUTOFF && ( noldaddconss != *naddconss || noldchgbds != *nchgbds ) )
3162  *result = SCIP_SUCCESS;
3163  }
3164 
3165  if ( *result != SCIP_CUTOFF && conshdlrdata->sdpconshdlrdata->diagzeroimplcuts )
3166  {
3167  noldaddconss = *naddconss;
3168  SCIP_CALL( diagZeroImpl(scip, conss, nconss, naddconss) );
3169  SCIPdebugMsg(scip, "Added %d cuts for implication from 0 diagonal.\n", *naddconss - noldaddconss);
3170  if ( noldaddconss != *naddconss )
3171  *result = SCIP_SUCCESS;
3172  }
3173 
3174  if ( *result != SCIP_CUTOFF && conshdlrdata->sdpconshdlrdata->twominorlinconss )
3175  {
3176  noldaddconss = *naddconss;
3177  SCIP_CALL( addTwoMinorLinConstraints(scip, conss, nconss, naddconss) );
3178  SCIPdebugMsg(scip, "Added %d linear constraints for 2 by 2 minors.\n", *naddconss - noldaddconss);
3179  if ( noldaddconss != *naddconss )
3180  *result = SCIP_SUCCESS;
3181  }
3182 
3183  if ( *result != SCIP_CUTOFF && conshdlrdata->sdpconshdlrdata->twominorprodconss )
3184  {
3185  noldaddconss = *naddconss;
3186  SCIP_CALL( addTwoMinorProdConstraints(scip, conss, nconss, naddconss) );
3187  SCIPdebugMsg(scip, "Added %d linear constraints for products of 2 by 2 minors.\n", *naddconss - noldaddconss);
3188  if ( noldaddconss != *naddconss )
3189  *result = SCIP_SUCCESS;
3190  }
3191  }
3192  }
3193 
3194  return SCIP_OKAY;
3195 }
3196 
3198 static
3199 SCIP_DECL_CONSTRANS(consTransSdp)
3200 {/*lint --e{715}*/
3201  SCIP_CONSDATA* sourcedata;
3202  SCIP_CONSDATA* targetdata;
3203 #ifdef OMP
3204  SCIP_CONSHDLRDATA* conshdlrdata;
3205 #endif
3206 #ifndef NDEBUG
3207  int snprintfreturn; /* used to check the return code of snprintf */
3208 #endif
3209  int i;
3210  char transname[SCIP_MAXSTRLEN];
3211 
3212  sourcedata = SCIPconsGetData(sourcecons);
3213  assert( sourcedata != NULL );
3214 
3215  SCIPdebugMsg(scip, "Transforming constraint <%s>\n", SCIPconsGetName(sourcecons));
3216 
3217 #ifdef OMP
3218  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3219  SCIPdebugMsg(scip, "Setting number of threads to %d via OpenMP in Openblas.\n", conshdlrdata->nthreads);
3220  omp_set_num_threads(conshdlrdata->nthreads);
3221 #endif
3222 
3223  SCIP_CALL( SCIPallocBlockMemory(scip, &targetdata) );
3224 
3225  /* copy some general data */
3226  targetdata->nvars = sourcedata->nvars;
3227  targetdata->nnonz = sourcedata->nnonz;
3228  targetdata->blocksize = sourcedata->blocksize;
3229 
3230  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(targetdata->nvarnonz), sourcedata->nvarnonz, sourcedata->nvars) );
3231 
3232  /* copy the non-constant nonzeros */
3233  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(targetdata->col), sourcedata->nvars) );
3234  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(targetdata->row), sourcedata->nvars) );
3235  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(targetdata->val), sourcedata->nvars) );
3236 
3237  for (i = 0; i < sourcedata->nvars; i++)
3238  {
3239  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(targetdata->col[i]), sourcedata->col[i], sourcedata->nvarnonz[i]) );
3240  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(targetdata->row[i]), sourcedata->row[i], sourcedata->nvarnonz[i]) );
3241  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(targetdata->val[i]), sourcedata->val[i], sourcedata->nvarnonz[i]) );
3242  }
3243  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(targetdata->vars), sourcedata->nvars) );
3244  if ( sourcedata->locks )
3245  {
3246  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(targetdata->locks), sourcedata->locks, sourcedata->nvars) );
3247  }
3248  else
3249  targetdata->locks = NULL;
3250 
3251  /* copy & transform the vars array */
3252  for (i = 0; i < sourcedata->nvars; i++)
3253  {
3254  targetdata->vars[i] = SCIPvarGetTransVar(sourcedata->vars[i]);
3255  SCIP_CALL( SCIPcaptureVar(scip, targetdata->vars[i]) );
3256  }
3257 
3258  /* copy the constant nonzeros */
3259  targetdata->constnnonz = sourcedata->constnnonz;
3260 
3261  if ( sourcedata->constnnonz > 0 )
3262  {
3263  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(targetdata->constcol), sourcedata->constcol, sourcedata->constnnonz));
3264  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(targetdata->constrow), sourcedata->constrow, sourcedata->constnnonz));
3265  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(targetdata->constval), sourcedata->constval, sourcedata->constnnonz));
3266  }
3267  else
3268  {
3269  targetdata->constcol = NULL;
3270  targetdata->constrow = NULL;
3271  targetdata->constval = NULL;
3272  }
3273 
3274  /* copy the maxrhsentry */
3275  targetdata->maxrhsentry = sourcedata->maxrhsentry;
3276 
3277  /* copy rankone */
3278  targetdata->rankone = sourcedata->rankone;
3279 
3280  /* copy maxevsubmat */
3281  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(targetdata->maxevsubmat), sourcedata->maxevsubmat, 2) );
3282 
3283  /* copy addedquadcons */
3284  targetdata->addedquadcons = sourcedata->addedquadcons;
3285 
3286  /* name the transformed constraint */
3287 #ifndef NDEBUG
3288  snprintfreturn = SCIPsnprintf(transname, SCIP_MAXSTRLEN, "t_%s", SCIPconsGetName(sourcecons));
3289  assert( snprintfreturn < SCIP_MAXSTRLEN ); /* check whether the name fits into the string */
3290 #else
3291  (void) SCIPsnprintf(transname, SCIP_MAXSTRLEN, "t_%s", SCIPconsGetName(sourcecons));
3292 #endif
3293 
3294  /* create target constraint */
3295  SCIP_CALL( SCIPcreateCons(scip, targetcons, transname, conshdlr, targetdata,
3296  SCIPconsIsInitial(sourcecons), SCIPconsIsSeparated(sourcecons), SCIPconsIsEnforced(sourcecons),
3297  SCIPconsIsChecked(sourcecons), SCIPconsIsPropagated(sourcecons), SCIPconsIsLocal(sourcecons),
3298  SCIPconsIsModifiable(sourcecons), SCIPconsIsDynamic(sourcecons), SCIPconsIsRemovable(sourcecons),
3299  SCIPconsIsStickingAtNode(sourcecons)) );
3300 
3301  return SCIP_OKAY;
3302 }
3303 
3305 static
3306 SCIP_DECL_CONSCHECK(consCheckSdp)
3307 {/*lint --e{715}*/
3308  SCIP_CONS* violcons;
3309  SCIP_CONSDATA* consdata;
3310  SCIP_CONSHDLRDATA* conshdlrdata;
3311  SCIP_VAR** linvars;
3312  SCIP_VAR* var;
3313  SCIP_SOL* bestrank1approx;
3314  SCIP_Real* matrix;
3315  SCIP_Real* fullmatrix;
3316  SCIP_Real* eigenvalues;
3317  SCIP_Real* eigenvectors;
3318  SCIP_Real* scaledeigenvectors;
3319  SCIP_Real* matrixC;
3320  SCIP_Real* matrixAj;
3321  SCIP_Real* linmatrix;
3322  SCIP_Real* rhsmatrix;
3323  SCIP_Real* lssolu;
3324  SCIP_Real* linvals;
3325  SCIP_Real* colmatrix;
3326  SCIP_Bool rank1result;
3327  SCIP_Bool stored;
3328 
3329  int c;
3330  int i;
3331  int j;
3332  int k;
3333  int l;
3334 #ifdef PRINTMATRICES
3335  int r;
3336  int s;
3337 #endif
3338  int idx;
3339  int blocksize;
3340  int nviolrank1 = 0;
3341  int nvars;
3342  int nrank1vars = 0;
3343  int linrows = 0; /* the number of rows of the linear equation system is given by the total number of entries
3344  (lower-triangular) in all violated rank-1 constraints altogether */
3345  int lincnt = 0;
3346  int nsdpvars;
3347  int* rank1considx;
3348  int* indviolrank1conss;
3349 
3350  SCIP_VAR** rank1consvars;
3351 
3352 
3353  assert( scip != NULL );
3354  assert( result != NULL );
3355  assert( conss != NULL );
3356 
3357  *result = SCIP_FEASIBLE;
3358 
3359 #ifdef PRINTMATRICES
3360  SCIP_CALL( SCIPprintSol(scip, sol, NULL, FALSE) );
3361 #endif
3362 
3363  /* check positive semidefiniteness */
3364  for (i = 0; i < nconss; ++i)
3365  {
3366  SCIP_CALL( SCIPconsSdpCheckSdpCons(scip, conss[i], sol, printreason, result) );
3367 #ifdef PRINTMATRICES
3368  SCIPinfoMessage(scip, NULL, "Solution is %d for constraint %s.\n", *result, SCIPconsGetName(conss[i]) );
3369 #endif
3370  if ( *result == SCIP_INFEASIBLE )
3371  return SCIP_OKAY;
3372  }
3373 
3374  /* check if heuristic should be executed */
3375  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3376  assert( conshdlrdata != NULL );
3377  if ( ! conshdlrdata->sdpconshdlrdata->rank1approxheur )
3378  {
3379  return SCIP_OKAY;
3380  }
3381 
3382  /* do not run in sub-SCIPs to avoid recursive reformulations due to rank 1 constraints */
3383  if ( SCIPgetSubscipDepth(scip) > 0 )
3384  return SCIP_OKAY;
3385 
3386  SCIP_CALL( SCIPallocBufferArray(scip, &indviolrank1conss, nconss) );
3387 
3388  /* check for rank-1 if necessary */
3389  SCIPdebugMsg(scip, "Check rank-1 constraints if there are any.\n");
3390  for (i = 0; i < nconss; ++i)
3391  {
3392  consdata = SCIPconsGetData(conss[i]);
3393  assert( consdata != NULL );
3394 
3395  if ( consdata->rankone )
3396  {
3397  SCIP_CALL( isMatrixRankOne(scip, conss[i], sol, &rank1result) );
3398  if ( ! rank1result )
3399  {
3400  /* save index of violated rank-1 constraint */
3401  indviolrank1conss[nviolrank1] = i;
3402  ++nviolrank1;
3403  }
3404  }
3405  }
3406 
3407  /* if there are no (violated) rank-1 constraint, we are finished. Otherwise, try to compute a feasible primal
3408  solution by computing the best rank-1 approximation for each violated rank-1 constraint and solve an LP to find a
3409  solution for the appearing variables */
3410  if ( nviolrank1 == 0 )
3411  {
3412  SCIPdebugMsg(scip, "Found no violated rank-1 constraints.\n");
3413  SCIPfreeBufferArray(scip, &indviolrank1conss);
3414  return SCIP_OKAY;
3415  }
3416 
3417  SCIPdebugMsg(scip, "Found %d violated rank-1 constraints, thus apply rank-1 approximation heuristic!\n", nviolrank1);
3418 
3419  /* we have to collect all variables appearing in violated SDPrank1-constraints first */
3420  nvars = SCIPgetNVars(scip);
3421  SCIP_CALL( SCIPallocBufferArray(scip, &rank1considx, nvars) );
3422  SCIP_CALL( SCIPallocBufferArray(scip, &rank1consvars, nvars) );
3423 
3424  for (j = 0; j < nvars; ++j)
3425  rank1considx[j] = -1;
3426 
3427  for (c = 0; c < nviolrank1; ++c)
3428  {
3429  assert( conss[indviolrank1conss[c]] != NULL );
3430 
3431  /* todo: write function to only get number of variables and variables of an SDP constraint */
3432  consdata = SCIPconsGetData(conss[indviolrank1conss[c]]);
3433  assert( consdata != NULL );
3434 
3435  nsdpvars = consdata->nvars;
3436  linrows += consdata->blocksize * (consdata->blocksize + 1) / 2;
3437 
3438  for (i = 0; i < nsdpvars; ++i)
3439  {
3440  var = consdata->vars[i];
3441  idx = SCIPvarGetProbindex(var);
3442  assert( 0 <= idx && idx < nvars );
3443  if ( rank1considx[idx] < 0 )
3444  {
3445  rank1consvars[nrank1vars] = var;
3446  rank1considx[idx] = nrank1vars++;
3447  }
3448  }
3449  }
3450 
3451  /* initialize matrix of linear equation system */
3452  SCIP_CALL( SCIPallocClearBufferArray(scip, &linmatrix, linrows * nrank1vars) );
3453  SCIP_CALL( SCIPallocClearBufferArray(scip, &rhsmatrix, MAX(linrows,nrank1vars)) );
3454 
3455  for (i = 0; i < nviolrank1; ++i)
3456  {
3457  /* get violated rank-1 constraint */
3458  violcons = conss[indviolrank1conss[i]];
3459  consdata = SCIPconsGetData(violcons);
3460 
3461  assert( consdata != NULL );
3462 
3463  blocksize = consdata->blocksize;
3464 
3465  SCIPdebugMsg(scip, "\n Start with violated rank-1 constraint %s, is %d out of %d violated rank-1 constraints.\n\n", SCIPconsGetName(violcons), i + 1, nviolrank1);
3466 
3467  /* allocate memory to store full matrix */
3468  SCIP_CALL( SCIPallocBufferArray(scip, &matrix, (blocksize * (blocksize+1))/2 ) );
3469  SCIP_CALL( SCIPallocBufferArray(scip, &fullmatrix, blocksize * blocksize ) );
3470  SCIP_CALL( SCIPallocBufferArray(scip, &eigenvalues, blocksize) );
3471  SCIP_CALL( SCIPallocBufferArray(scip, &eigenvectors, blocksize * blocksize) );
3472  SCIP_CALL( SCIPallocBufferArray(scip, &matrixC, blocksize * blocksize) );
3473  SCIP_CALL( SCIPallocBufferArray(scip, &matrixAj, blocksize * blocksize) );
3474  SCIP_CALL( SCIPallocBufferArray(scip, &linvars, consdata->nvars) );
3475  SCIP_CALL( SCIPallocBufferArray(scip, &linvals, consdata->nvars) );
3476 
3477  /* compute the matrix \f$ \sum_j A_j y_j - A_0 \f$ */
3478  SCIP_CALL( computeSdpMatrix(scip, violcons, sol, matrix) );
3479 
3480  /* expand it because LAPACK wants the full matrix instead of the lower triangular part */
3481  SCIP_CALL( expandSymMatrix(blocksize, matrix, fullmatrix) );
3482 
3483 #ifdef PRINTMATRICES
3484  /* SCIPSDP uses row-first format! */
3485  printf("Full SDP-constraint matrix Z: \n");
3486  for (j = 0; j < blocksize; ++j)
3487  {
3488  for (k = 0; k < blocksize; ++k)
3489  printf("%.5f ", fullmatrix[j*blocksize + k]);
3490  printf("\n");
3491  }
3492 
3493  /* Double-check that fullmatrix is in row-first format */
3494  printf("Full SDP-constraint matrix Z in row-first format: \n");
3495  for (j = 0; j < blocksize * blocksize; ++j)
3496  printf("%.5f ", fullmatrix[j]);
3497  printf("\n");
3498 #endif
3499 
3500  /* compute EVD */
3501  SCIP_CALL( SCIPlapackComputeEigenvectorDecomposition(SCIPbuffer(scip), blocksize, fullmatrix, eigenvalues, eigenvectors) );
3502 
3503 #ifdef PRINTMATRICES
3504  /* caution: LAPACK uses column-first format! */
3505  printf("Eigenvectors of Z: \n");
3506  for (j = 0; j < blocksize; ++j)
3507  {
3508  for (k = 0; k < blocksize; ++k)
3509  printf("%.5f ", eigenvectors[k*blocksize + j]);
3510  printf("\n");
3511  }
3512 
3513  /* Double-check that eigenvectors is in column-first format */
3514  printf("Eigenvectors of Z in column-first format: \n");
3515  for (j = 0; j < blocksize * blocksize; ++j)
3516  printf("%.5f ", eigenvectors[j]);
3517  printf("\n");
3518 
3519  printf("Eigenvalues of Z: \n");
3520  for (j = 0; j < blocksize; ++j)
3521  printf("%.5f ", eigenvalues[j]);
3522  printf("\n");
3523 #endif
3524 
3525  /* duplicate memory of eigenvectors to compute diag(0,...,0,lambda_max) * U^T */
3526  SCIP_CALL( SCIPduplicateBufferArray(scip, &scaledeigenvectors, eigenvectors, blocksize*blocksize) );
3527 
3528  /* set all eigenvalues except the largest to zero (using the property that LAPACK returns them in ascending
3529  order) */
3530  for (j = 0; j < blocksize-1; ++j)
3531  {
3532  assert( ! SCIPisFeasNegative(scip, eigenvalues[j]) ); /* otherwise constraint is not psd */
3533  eigenvalues[j] = 0.0;
3534  }
3535 
3536  /* compute diag(0,...,0,lambda_max) * U^T. Note that scaledeigenvectors on entry is U in column-first format,
3537  i.e., U^T in row-first format. Thus, on exit, scaledeigenvectors contains diag(0,...,0,lambda_max) * U^T in
3538  row-first format. */
3539  SCIP_CALL( scaleRowsMatrix(blocksize, scaledeigenvectors, eigenvalues) );
3540 
3541 #ifdef PRINTMATRICES
3542  printf("Scaled eigenvectors of Z (only keep largest eigenvalue and corresponding eigenvector) : \n");
3543  for (j = 0; j < blocksize; ++j)
3544  {
3545  for (k = 0; k < blocksize; ++k)
3546  {
3547  printf("%.5f ", scaledeigenvectors[j*blocksize + k]);
3548  }
3549  printf("\n");
3550  }
3551 
3552  /* Double-check that scaledeigenvectors is in row-first format */
3553  printf("Scaled eigenvectors of Z in row-first format: \n");
3554  for (j = 0; j < blocksize * blocksize; ++j)
3555  printf("%.5f ", scaledeigenvectors[j]);
3556  printf("\n");
3557 #endif
3558 
3559  /* compute U * [diag(0,...,0,lambda_max) * U^T]: Since eigenvectors, which contains U, already comes in LAPACK's
3560  column-first format, eigenvectors does not need to be transposed! scaledeigenvectors models
3561  [diag(0,...,0,lambda_max) * U^T in SCIPSDP's row-first format, so that scaledeigenvectors needs to be
3562  transposed for LAPACK! */
3563  SCIP_CALL( SCIPlapackMatrixMatrixMult(blocksize, blocksize, eigenvectors, FALSE, blocksize, blocksize, scaledeigenvectors,
3564  TRUE, fullmatrix) );
3565 
3566 #ifdef PRINTMATRICES
3567  printf("Best rank-1 approximation of Z: \n");
3568  for (j = 0; j < blocksize; ++j)
3569  {
3570  for (k = 0; k < blocksize; ++k)
3571  printf("%.5f ", fullmatrix[j*blocksize + k]);
3572  printf("\n");
3573  }
3574 
3575  /* Double-check that fullmatrix (best rank-1 approximation) is in row-first format */
3576  printf("Best rank-1 approximation of Z in row-first format: \n");
3577  for (j = 0; j < blocksize * blocksize; ++j)
3578  printf("%.5f ", fullmatrix[j]);
3579  printf("\n");
3580 #endif
3581 
3582  /* update linear equation system */
3583 
3584  /* compute constant matrix A_0 in row-first format*/
3585  SCIP_CALL( SCIPconsSdpGetFullConstMatrix(scip, violcons, matrixC) );
3586 
3587 #ifdef PRINTMATRICES
3588  printf("Constant matrix A_0 of SDP-constraint: \n");
3589  for (j = 0; j < blocksize; ++j)
3590  {
3591  for (k = 0; k < blocksize; ++k)
3592  printf("%.5f ", matrixC[j*blocksize + k]);
3593  printf("\n");
3594  }
3595 
3596  /* Double-check that matrixA0 is in row-first format */
3597  printf("Constant matrix A_0 of SDP-constraint in row-first format: \n");
3598  for (j = 0; j < blocksize * blocksize; ++j)
3599  printf("%.5f ", matrixC[j]);
3600  printf("\n");
3601 #endif
3602 
3603  for (j = 0; j < blocksize; ++j)
3604  {
3605  for (k = 0; k <= j; ++k)
3606  {
3607  for (l = 0; l < consdata->nvars; ++l)
3608  {
3609  /* compute matrix A_j in row-first format */
3610  SCIP_CALL( SCIPconsSdpGetFullAj(scip, violcons, l, matrixAj) );
3611 
3612 #ifdef PRINTMATRICES
3613  printf("Coefficient matrix A_%d of SDP-constraint: \n", l+1);
3614  for (r = 0; r < blocksize; ++r)
3615  {
3616  for (s = 0; s < blocksize; ++s)
3617  printf("%.5f ", matrixAj[r*blocksize + s]);
3618  printf("\n");
3619  }
3620 
3621  /* Double-check that matrixAj is in row-first format */
3622  printf("Constant matrix A_0 of SDP-constraint in row-first format: \n");
3623  for (r = 0; r < blocksize * blocksize; ++r)
3624  printf("%.5f ", matrixAj[r]);
3625  printf("\n");
3626 #endif
3627 
3628  idx = SCIPvarGetProbindex(consdata->vars[l]);
3629  idx = rank1considx[idx];
3630  assert( 0 <= idx && idx < nrank1vars );
3631  assert( lincnt <= linrows );
3632  linmatrix[lincnt * nrank1vars + idx] = matrixAj[j * blocksize + k];
3633  }
3634  rhsmatrix[lincnt] = matrixC[j * blocksize + k] + fullmatrix[j * blocksize + k];
3635  ++lincnt;
3636  }
3637  }
3638 
3639  /* free memory for full matrix, eigenvalues and eigenvectors */
3640  SCIPfreeBufferArray(scip, &linvals);
3641  SCIPfreeBufferArray(scip, &linvars);
3642  SCIPfreeBufferArray(scip, &matrixAj);
3643  SCIPfreeBufferArray(scip, &matrixC);
3644  SCIPfreeBufferArray(scip, &scaledeigenvectors);
3645  SCIPfreeBufferArray(scip, &eigenvectors);
3646  SCIPfreeBufferArray(scip, &eigenvalues);
3647  SCIPfreeBufferArray(scip, &fullmatrix);
3648  SCIPfreeBufferArray(scip, &matrix);
3649  }
3650 
3651  assert( lincnt == linrows );
3652 
3653 #ifdef PRINTMATRICES
3654  printf("Matrix for linear equation system, in row-first format:\n");
3655  for (j = 0; j < linrows; ++j)
3656  {
3657  for (k = 0; k < nrank1vars; ++k)
3658  {
3659  printf("%.5f ", linmatrix[j * nrank1vars + k]);
3660  }
3661  printf("\n");
3662  }
3663 
3664  /* Double-check that linmatrix is in row-first format */
3665  printf("Matrix for linear equation system in row-first format: \n");
3666  for (r = 0; r < linrows * nrank1vars; ++r)
3667  printf("%.5f ", linmatrix[r]);
3668  printf("\n");
3669 
3670  printf("Right-hand for linear equation system:\n");
3671  for (j = 0; j < nrank1vars; ++j)
3672  {
3673  printf("%.5f ", rhsmatrix[j]);
3674  }
3675  printf("\n");
3676 #endif
3677 
3678  /* solve linear equation system with LAPACK */
3679  SCIP_CALL( SCIPallocBufferArray(scip, &lssolu, nrank1vars) );
3680 
3681  /* caution: LAPACK wants matrices in columns-first format, but SCIPSDP represents matrices in row-first format */
3682  SCIP_CALL( SCIPallocBufferArray(scip, &colmatrix, linrows * nrank1vars ) );
3683 
3684  SCIP_CALL( convertRowToColFormatFullMatrix(linrows, nrank1vars, linmatrix, colmatrix) );
3685 
3686 #ifdef PRINTMATRICES
3687  printf("Matrix for linear equation system, in col-first format:\n");
3688  for (j = 0; j < linrows; ++j)
3689  {
3690  for (l = 0; l < nrank1vars; ++l)
3691  {
3692  printf("%.5f ", colmatrix[l * linrows + j]);
3693  }
3694  printf("\n");
3695  }
3696 
3697  /* Double-check that colmatrix is in col-first format */
3698  printf("Matrix for linear equation system in col-first format: \n");
3699  for (r = 0; r < linrows * nrank1vars; ++r)
3700  printf("%.5f ", colmatrix[r]);
3701  printf("\n");
3702 #endif
3703 
3704  SCIP_CALL( SCIPlapackLinearSolve( SCIPbuffer(scip), linrows, nrank1vars, colmatrix, rhsmatrix, lssolu) );
3705 
3706  /* copy current solution */
3707  SCIP_CALL( SCIPcreateSolCopy(scip, &bestrank1approx, sol) );
3708 
3709  /* update solution with values from linear equations system solution from LAPACK */
3710  for (i = 0; i < nrank1vars; ++i)
3711  {
3712  var = rank1consvars[i];
3713  SCIP_CALL( SCIPsetSolVal(scip, bestrank1approx, var, lssolu[i]) );
3714  }
3715 
3716  SCIP_CALL( SCIPtrySolFree(scip, &bestrank1approx, FALSE, TRUE, TRUE, TRUE, TRUE, &stored) );
3717  if ( stored )
3718  SCIPdebugMsg(scip, "Best Rank-1 Approximation Heuristic found feasible primal solution\n");
3719  else
3720  SCIPdebugMsg(scip, "Primal solution found by Best Rank-1 Approximation Heuristic is not feasible!\n");
3721 
3722  SCIPfreeBufferArray(scip, &colmatrix);
3723  SCIPfreeBufferArray(scip, &lssolu);
3724  SCIPfreeBufferArray(scip, &rhsmatrix);
3725  SCIPfreeBufferArray(scip, &linmatrix);
3726  SCIPfreeBufferArray(scip, &rank1consvars);
3727  SCIPfreeBufferArray(scip, &rank1considx);
3728  SCIPfreeBufferArray(scip, &indviolrank1conss);
3729 
3730  return SCIP_OKAY;
3731 }
3732 
3737 static
3738 SCIP_DECL_CONSENFOPS(consEnfopsSdp)
3739 {/*lint --e{715}*/
3740  int i;
3741 
3742  assert( scip != NULL );
3743  assert( result != NULL );
3744  assert( conss != NULL );
3745 
3746  *result = SCIP_DIDNOTRUN;
3747 
3748  if ( objinfeasible )
3749  {
3750  SCIPdebugMsg(scip, "-> pseudo solution is objective infeasible, return.\n");
3751  return SCIP_OKAY;
3752  }
3753 
3754  for (i = 0; i < nconss; ++i)
3755  {
3756  SCIP_CALL( SCIPconsSdpCheckSdpCons(scip, conss[i], NULL, FALSE, result) );
3757 
3758  if (*result == SCIP_INFEASIBLE)
3759  {
3760  /* if it is infeasible for one SDP constraint, it is infeasible for the whole problem */
3761  SCIPdebugMsg(scip, "-> pseudo solution infeasible for SDP-constraint %s, return.\n", SCIPconsGetName(conss[i]));
3762  return SCIP_OKAY;
3763  }
3764  }
3765 
3766  *result = SCIP_FEASIBLE;
3767 
3768  SCIPdebugMsg(scip, "-> pseudo solution feasible for all SDP-constraints.\n");
3769 
3770  return SCIP_OKAY;
3771 }
3772 
3773 
3776 static
3777 SCIP_DECL_CONSENFOLP(consEnfolpSdp)
3778 {/*lint --e{715}*/
3779  assert( scip != NULL );
3780  assert( conshdlr != NULL );
3781  assert( conss != NULL );
3782  assert( result != NULL );
3783 
3784  return enforceConstraint(scip, conshdlr, conss, nconss, NULL, result);
3785 }
3786 
3789 static
3790 SCIP_DECL_CONSENFORELAX(consEnforelaxSdp)
3791 {/*lint --e{715}*/
3792  assert( scip != NULL );
3793  assert( conshdlr != NULL );
3794  assert( conss != NULL );
3795  assert( result != NULL );
3796 
3797  return enforceConstraint(scip, conshdlr, conss, nconss, sol, result);
3798 }
3799 
3801 static
3802 SCIP_DECL_CONSSEPASOL(consSepasolSdp)
3803 {/*lint --e{715}*/
3804  int i;
3805 
3806  assert( result != NULL );
3807  *result = SCIP_DIDNOTFIND;
3808 
3809  for (i = 0; i < nusefulconss; ++i)
3810  {
3811  SCIP_CALL( separateSol(scip, conshdlr, conss[i], sol, result) );
3812  }
3813 
3814  return SCIP_OKAY;
3815 }
3816 
3818 static
3819 SCIP_DECL_CONSSEPALP(consSepalpSdp)
3820 {/*lint --e{715}*/
3821  int i;
3822 
3823  assert( result != NULL );
3824  *result = SCIP_DIDNOTFIND;
3825 
3826  for (i = 0; i < nusefulconss; ++i)
3827  {
3828  SCIP_CALL( separateSol(scip, conshdlr, conss[i], NULL, result) );
3829  }
3830 
3831  return SCIP_OKAY;
3832 }
3833 
3835 static
3836 SCIP_DECL_CONSDELETE(consDeleteSdp)
3837 {/*lint --e{715}*/
3838  int i;
3839 
3840  assert( cons != NULL );
3841  assert( consdata != NULL );
3842 
3843  SCIPdebugMsg(scip, "deleting SDP constraint <%s>.\n", SCIPconsGetName(cons));
3844 
3845  /* release memory for rank one constraint */
3846  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->maxevsubmat, 2);
3847 
3848  for (i = 0; i < (*consdata)->nvars; i++)
3849  {
3850  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->val[i], (*consdata)->nvarnonz[i]);
3851  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->row[i], (*consdata)->nvarnonz[i]);
3852  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->col[i], (*consdata)->nvarnonz[i]);
3853  }
3854 
3855  /* release all variables */
3856  for (i = 0; i < (*consdata)->nvars; i++)
3857  {
3858  SCIP_CALL( SCIPreleaseVar(scip, &((*consdata)->vars[i])) );
3859  }
3860 
3861  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->vars, (*consdata)->nvars);
3862  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->locks, (*consdata)->nvars);
3863  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->constval, (*consdata)->constnnonz);
3864  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->constrow, (*consdata)->constnnonz);
3865  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->constcol, (*consdata)->constnnonz);
3866  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->val, (*consdata)->nvars);
3867  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->row, (*consdata)->nvars);
3868  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->col, (*consdata)->nvars);
3869  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->nvarnonz, (*consdata)->nvars);
3870  SCIPfreeBlockMemory(scip, consdata);
3871 
3872  return SCIP_OKAY;
3873 }
3874 
3876 static
3877 SCIP_DECL_CONSFREE(consFreeSdp)
3878 {/*lint --e{715}*/
3879  SCIP_CONSHDLRDATA* conshdlrdata;
3880 
3881  SCIPdebugMsg(scip, "Freeing constraint handler <%s>.\n", SCIPconshdlrGetName(conshdlr));
3882 
3883  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3884  assert( conshdlrdata != NULL );
3885 
3886  SCIPfreeMemory(scip, &conshdlrdata);
3887  SCIPconshdlrSetData(conshdlr, NULL);
3888 
3889  return SCIP_OKAY;
3890 }
3891 
3893 static
3894 SCIP_DECL_CONSHDLRCOPY(conshdlrCopySdp)
3895 {
3896  assert(scip != NULL);
3897  assert(conshdlr != NULL);
3898  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
3899 
3900  SCIP_CALL( SCIPincludeConshdlrSdp(scip) );
3901 
3902  *valid = TRUE;
3903 
3904  return SCIP_OKAY;
3905 }
3906 
3908 static
3909 SCIP_DECL_CONSHDLRCOPY(conshdlrCopySdpRank1)
3910 {
3911  assert(scip != NULL);
3912  assert(conshdlr != NULL);
3913  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLRRANK1_NAME) == 0);
3914 
3915  SCIP_CALL( SCIPincludeConshdlrSdpRank1(scip) );
3916 
3917  *valid = TRUE;
3918 
3919  return SCIP_OKAY;
3920 }
3921 
3923 static
3924 SCIP_DECL_CONSCOPY(consCopySdp)
3925 {/*lint --e{715}*/
3926  char copyname[SCIP_MAXSTRLEN];
3927  SCIP_CONSDATA* sourcedata;
3928  SCIP_Bool success;
3929  SCIP_VAR** targetvars;
3930  SCIP_VAR* var;
3931  int i;
3932 #ifndef NDEBUG
3933  int snprintfreturn; /* used to check the return code of snprintf */
3934 #endif
3935 
3936  assert( scip != NULL );
3937  assert( sourcescip != NULL );
3938  assert( sourcecons != NULL );
3939  assert( valid != NULL );
3940 
3941  SCIPdebugMsg(scip, "Copying SDP constraint <%s>\n", SCIPconsGetName(sourcecons));
3942 
3943  *valid = TRUE;
3944 
3945  /* as we can only map active variables, we have to make sure, that the constraint contains no fixed or (multi-)aggregated vars, after
3946  * exitpresolve (stage 6) this should always be the case, earlier than that we need to call fixAndAggrVars */
3947  if ( SCIPgetStage(sourcescip) <= SCIP_STAGE_EXITPRESOLVE )
3948  {
3949  SCIP_CALL( fixAndAggrVars(sourcescip, &sourcecons, 1, TRUE) );
3950  }
3951 
3952  sourcedata = SCIPconsGetData(sourcecons);
3953  assert( sourcedata != NULL );
3954  assert( sourcedata->rankone || strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(sourcecons)), CONSHDLR_NAME) == 0 );
3955  assert( ! sourcedata->rankone || strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(sourcecons)), CONSHDLRRANK1_NAME) == 0 );
3956 
3957  SCIP_CALL( SCIPallocBufferArray(scip, &targetvars, sourcedata->nvars) );
3958 
3959  /* map all variables in the constraint */
3960  for (i = 0; i < sourcedata->nvars; i++)
3961  {
3962  SCIP_CALL( SCIPgetVarCopy(sourcescip, scip, sourcedata->vars[i], &var, varmap, consmap, global, &success) );
3963  if ( success )
3964  targetvars[i] = var;
3965  else
3966  *valid = FALSE;
3967  }
3968 
3969  /* name the copied constraint */
3970 #ifndef NDEBUG
3971  snprintfreturn = SCIPsnprintf(copyname, SCIP_MAXSTRLEN, "c_%s", name == NULL ? SCIPconsGetName(sourcecons) : name);
3972  assert( snprintfreturn < SCIP_MAXSTRLEN ); /* check whether the name fits into the string */
3973 #else
3974  (void) SCIPsnprintf(copyname, SCIP_MAXSTRLEN, "c_%s", name == NULL ? SCIPconsGetName(sourcecons) : name);
3975 #endif
3976 
3977  /* create the new constraint */
3978  if ( ! sourcedata->rankone )
3979  {
3980  SCIP_CALL( SCIPcreateConsSdp(scip, cons, copyname, sourcedata->nvars, sourcedata->nnonz, sourcedata->blocksize, sourcedata->nvarnonz,
3981  sourcedata->col, sourcedata->row, sourcedata->val, targetvars, sourcedata->constnnonz,
3982  sourcedata->constcol, sourcedata->constrow, sourcedata->constval) );
3983  }
3984  else
3985  {
3986  SCIP_CALL( SCIPcreateConsSdpRank1(scip, cons, copyname, sourcedata->nvars, sourcedata->nnonz, sourcedata->blocksize, sourcedata->nvarnonz,
3987  sourcedata->col, sourcedata->row, sourcedata->val, targetvars, sourcedata->constnnonz,
3988  sourcedata->constcol, sourcedata->constrow, sourcedata->constval) );
3989  }
3990 
3991  SCIPfreeBufferArray(scip, &targetvars);
3992 
3993  return SCIP_OKAY;
3994 }
3995 
3997 static
3998 SCIP_DECL_CONSPRINT(consPrintSdp)
3999 {/*lint --e{715}*/
4000 #ifdef PRINT_HUMAN_READABLE
4001  SCIP_CONSDATA* consdata;
4002  SCIP_Real* fullmatrix;
4003  int v;
4004  int i;
4005  int j;
4006 
4007  assert( scip != NULL );
4008  assert( cons != NULL );
4009 
4010  consdata = SCIPconsGetData(cons);
4011  assert( consdata != NULL );
4012 
4013  SCIP_CALL( SCIPallocBufferArray(scip, &fullmatrix, consdata->blocksize * consdata->blocksize) );
4014 
4015  /* print rank1 information */
4016  SCIPinfoMessage(scip, file, "rank-1? %d\n", consdata->rankone);
4017 
4018  /* print the non-constant matrices, for this they first have to be assembled in fullmatrix */
4019  for (v = 0; v < consdata->nvars; v++)
4020  {
4021  /* assemble the matrix */
4022 
4023  /* first initialize it with zero */
4024  for (i = 0; i < consdata->blocksize; i++)
4025  {
4026  for (j = 0; j < consdata->blocksize; j++)
4027  fullmatrix[i * consdata->blocksize + j] = 0.0; /*lint !e679*/
4028  }
4029 
4030  /* then add the nonzeros */
4031  for (i = 0; i < consdata->nvarnonz[v]; i++)
4032  {
4033  fullmatrix[consdata->row[v][i] * consdata->blocksize + consdata->col[v][i]] = consdata->val[v][i]; /* lower triangular entry */ /*lint !e679*/
4034  fullmatrix[consdata->col[v][i] * consdata->blocksize + consdata->row[v][i]] = consdata->val[v][i]; /* upper triangular entry */ /*lint !e679*/
4035  }
4036 
4037  /* print it */
4038  SCIPinfoMessage(scip, file, "+\n");
4039  for (i = 0; i < consdata->blocksize; i++)
4040  {
4041  SCIPinfoMessage(scip, file, "( ");
4042  for (j = 0; j < consdata->blocksize; j++)
4043  SCIPinfoMessage(scip, file, "%g ", fullmatrix[i * consdata->blocksize + j]); /*lint !e679*/
4044  SCIPinfoMessage(scip, file, ")\n");
4045  }
4046  SCIPinfoMessage(scip, file, "* %s\n", SCIPvarGetName(consdata->vars[v]));
4047  }
4048 
4049  /* print the constant-matrix */
4050 
4051  /* assemble the matrix */
4052 
4053  /* first initialize it with zero */
4054  for (i = 0; i < consdata->blocksize; i++)
4055  {
4056  for (j = 0; j < consdata->blocksize; j++)
4057  fullmatrix[i * consdata->blocksize + j] = 0.0; /*lint !e679*/
4058  }
4059 
4060  /* then add the nonzeros */
4061  for (i = 0; i < consdata->constnnonz; i++)
4062  {
4063  fullmatrix[consdata->constrow[i] * consdata->blocksize + consdata->constcol[i]] = consdata->constval[i]; /* lower triangular entry */ /*lint !e679*/
4064  fullmatrix[consdata->constcol[i] * consdata->blocksize + consdata->constrow[i]] = consdata->constval[i]; /* upper triangular entry */ /*lint !e679*/
4065  }
4066 
4067  /* print it */
4068  SCIPinfoMessage(scip, file, "-\n");
4069  for (i = 0; i < consdata->blocksize; i++)
4070  {
4071  SCIPinfoMessage(scip, file, "( ");
4072  for (j = 0; j < consdata->blocksize; j++)
4073  SCIPinfoMessage(scip, file, "%g ", fullmatrix[i * consdata->blocksize + j]); /*lint !e679*/
4074  SCIPinfoMessage(scip, file, ")\n");
4075  }
4076  SCIPinfoMessage(scip, file, ">= 0\n");
4077 
4078  /* print rank1 information */
4079  SCIPinfoMessage(scip, file, "rank-1? %d\n", consdata->rankone);
4080 
4081  SCIPfreeBufferArray(scip, &fullmatrix);
4082 
4083  return SCIP_OKAY;
4084 #else
4085  SCIP_CONSDATA* consdata;
4086  int i;
4087  int v;
4088 
4089  assert( scip != NULL );
4090  assert( cons != NULL );
4091 
4092  consdata = SCIPconsGetData(cons);
4093 
4094  /* print blocksize */
4095  SCIPinfoMessage(scip, file, "%d\n", consdata->blocksize);
4096 
4097  /* print rank1 information */
4098  SCIPinfoMessage(scip, file, " rank-1? %d\n", consdata->rankone);
4099 
4100  /* print A_0 if it exists */
4101  if ( consdata->constnnonz > 0 )
4102  {
4103  SCIPinfoMessage(scip, file, " A_0: ");
4104 
4105  for (i = 0; i < consdata->constnnonz; i++)
4106  {
4107  if ( i < consdata->constnnonz - 1 )
4108  SCIPinfoMessage(scip, file, "(%d,%d):%.15g, ", consdata->constrow[i], consdata->constcol[i], consdata->constval[i]);
4109  else
4110  SCIPinfoMessage(scip, file, "(%d,%d):%.15g", consdata->constrow[i], consdata->constcol[i], consdata->constval[i]);
4111  }
4112  SCIPinfoMessage(scip, file, "\n");
4113  }
4114 
4115  /* print other matrices */
4116  for (v = 0; v < consdata->nvars; v++)
4117  {
4118  SCIPinfoMessage(scip, file, " <%s>: ", SCIPvarGetName(consdata->vars[v]));
4119  for (i = 0; i < consdata->nvarnonz[v]; i++)
4120  {
4121  if ( i < consdata->nvarnonz[v] - 1 || v < consdata->nvars - 1 )
4122  SCIPinfoMessage(scip, file, "(%d,%d):%.15g, ", consdata->row[v][i], consdata->col[v][i], consdata->val[v][i]);
4123  else
4124  SCIPinfoMessage(scip, file, "(%d,%d):%.15g", consdata->row[v][i], consdata->col[v][i], consdata->val[v][i]);
4125  }
4126  /* if this is not the last variable, add a newline */
4127  if (v < consdata->nvars - 1)
4128  {
4129  SCIPinfoMessage(scip, file, "\n");
4130  }
4131  }
4132 
4133  return SCIP_OKAY;
4134 #endif
4135 }
4136 
4138 static
4139 SCIP_DECL_CONSPARSE(consParseSdp)
4140 { /*lint --e{715}*/
4141  SCIP_Bool parsesuccess;
4142  SCIP_CONSDATA* consdata = NULL;
4143  char* pos;
4144  int currentsize;
4145  int nvars;
4146  int i;
4147  int v;
4148  int rankoneint;
4149 
4150  assert( scip != NULL );
4151  assert( str != NULL );
4152 
4153  nvars = SCIPgetNVars(scip);
4154 
4155  assert( success != NULL );
4156  *success = TRUE;
4157 
4158  /* create constraint data */
4159  SCIP_CALL( SCIPallocBlockMemory(scip, &consdata) );
4160  consdata->nvars = 0;
4161  consdata->nnonz = 0;
4162  consdata->constnnonz = 0;
4163  consdata->rankone = 0;
4164  consdata->addedquadcons = FALSE;
4165  consdata->locks = NULL;
4166 
4167  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->nvarnonz, nvars) );
4168  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->col, nvars) );
4169  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->row, nvars) );
4170  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->val, nvars) );
4171  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->vars, nvars));
4172 
4173  consdata->constcol = NULL;
4174  consdata->constrow = NULL;
4175  consdata->constval = NULL;
4176 
4177  /* parse the blocksize */
4178  parsesuccess = SCIPstrToIntValue(str, &(consdata->blocksize), &pos);
4179  *success = *success && parsesuccess;
4180 
4181  /* skip whitespace */
4182  while ( isspace((unsigned char)*pos) )
4183  pos++;
4184 
4185  /* parse the rank1-information */
4186  if ( pos[0] == 'r' && pos[1] == 'a' && pos[2] == 'n' && pos[3] == 'k' && pos[4] == '-' && pos[5] == '1' && pos[6] == '?' )
4187  {
4188  pos += 8; /* we skip "rank1-? " */
4189  parsesuccess = SCIPstrToIntValue(pos, &rankoneint, &pos);
4190  consdata->rankone = (SCIP_Bool) rankoneint;
4191  *success = *success && parsesuccess;
4192 
4193  printf("rank-1? %d\n", rankoneint);
4194  }
4195 
4196  /* skip whitespace */
4197  while( isspace((unsigned char)*pos) )
4198  pos++;
4199 
4200  /* check if there is a constant part */
4201  if ( pos[0] == 'A' && pos[1] == '_' && pos[2] == '0' )
4202  {
4203  pos += 5; /* we skip "A_0: " */
4204 
4205  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->constcol, PARSE_STARTSIZE) );
4206  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->constrow, PARSE_STARTSIZE) );
4207  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->constval, PARSE_STARTSIZE) );
4208 
4209  currentsize = PARSE_STARTSIZE;
4210 
4211  /* as long as there is another entry for the constant part, parse it */
4212  while (pos[0] == '(')
4213  {
4214  pos++; /* remove the '(' */
4215 
4216  /* check if we need to enlarge the arrays */
4217  if ( consdata->constnnonz == currentsize )
4218  {
4219  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->constcol, currentsize, PARSE_SIZEFACTOR * currentsize) );
4220  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->constrow, currentsize, PARSE_SIZEFACTOR * currentsize) );
4221  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->constval, currentsize, PARSE_SIZEFACTOR * currentsize) );
4222  currentsize *= PARSE_SIZEFACTOR;
4223  }
4224 
4225  parsesuccess = SCIPstrToIntValue(pos, &(consdata->constrow[consdata->constnnonz]), &pos);
4226  *success = *success && parsesuccess;
4227  assert( consdata->constrow[consdata->constnnonz] < consdata->blocksize );
4228  pos++; /* remove the ',' */
4229  parsesuccess = SCIPstrToIntValue(pos, &(consdata->constcol[consdata->constnnonz]), &pos);
4230  *success = *success && parsesuccess;
4231  assert( consdata->constcol[consdata->constnnonz] < consdata->blocksize );
4232  pos += 2; /* remove the "):" */
4233  parsesuccess = SCIPstrToRealValue(pos, &(consdata->constval[consdata->constnnonz]), &pos);
4234  *success = *success && parsesuccess;
4235  pos ++; /* remove the "," */
4236 
4237  /* if we got an entry in the upper triangular part, switch the entries for lower triangular */
4238  if ( consdata->constcol[consdata->constnnonz] > consdata->constrow[consdata->constnnonz] )
4239  {
4240  i = consdata->constcol[consdata->constnnonz];
4241  consdata->constcol[consdata->constnnonz] = consdata->constrow[consdata->constnnonz];
4242  consdata->constrow[consdata->constnnonz] = i;
4243  }
4244 
4245  consdata->constnnonz++;
4246 
4247  /* skip whitespace */
4248  while( isspace((unsigned char)*pos) )
4249  pos++;
4250  }
4251 
4252  /* resize the arrays to their final size */
4253  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->constcol, currentsize, consdata->constnnonz) );
4254  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->constrow, currentsize, consdata->constnnonz) );
4255  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->constval, currentsize, consdata->constnnonz) );
4256  }
4257 
4258  /* skip whitespace */
4259  while( isspace((unsigned char)*pos) )
4260  pos++;
4261 
4262  /* parse the non-constant part */
4263 
4264  /* while there is another variable, parse it */
4265  while ( pos[0] == '<' )
4266  {
4267  /* add the variable to consdata->vars and create the corresponding nonzero arrays */
4268  SCIP_CALL( SCIPparseVarName(scip, pos, &(consdata->vars[consdata->nvars]), &pos) );
4269  SCIP_CALL( SCIPcaptureVar(scip, consdata->vars[consdata->nvars]) );
4270 
4271  consdata->nvarnonz[consdata->nvars] = 0;
4272  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(consdata->col[consdata->nvars]), PARSE_STARTSIZE));
4273  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(consdata->row[consdata->nvars]), PARSE_STARTSIZE));
4274  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(consdata->val[consdata->nvars]), PARSE_STARTSIZE));
4275  consdata->nvars++;
4276  currentsize = PARSE_STARTSIZE;
4277 
4278  pos += 2; /* remove the ": " */
4279 
4280  /* while there is another entry, parse it */
4281  while (pos[0] == '(')
4282  {
4283  pos++; /* remove the '(' */
4284 
4285  /* check if we need to enlarge the arrays */
4286  if ( consdata->nvarnonz[consdata->nvars - 1] == currentsize )
4287  {
4288  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->col[consdata->nvars - 1], currentsize, PARSE_SIZEFACTOR * currentsize) );
4289  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->row[consdata->nvars - 1], currentsize, PARSE_SIZEFACTOR * currentsize) );
4290  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->val[consdata->nvars - 1], currentsize, PARSE_SIZEFACTOR * currentsize) );
4291  currentsize *= PARSE_SIZEFACTOR;
4292  }
4293 
4294  parsesuccess = SCIPstrToIntValue(pos, &(consdata->row[consdata->nvars - 1][consdata->nvarnonz[consdata->nvars - 1]]), &pos);
4295  *success = *success && parsesuccess;
4296  assert( consdata->row[consdata->nvars - 1][consdata->nvarnonz[consdata->nvars - 1]] < consdata->blocksize );
4297  pos++; /* remove the ',' */
4298  parsesuccess = SCIPstrToIntValue(pos, &(consdata->col[consdata->nvars - 1][consdata->nvarnonz[consdata->nvars - 1]]), &pos);
4299  *success = *success && parsesuccess;
4300  assert( consdata->col[consdata->nvars - 1][consdata->nvarnonz[consdata->nvars - 1]] < consdata->blocksize );
4301  pos += 2; /* remove the "):" */
4302  parsesuccess = SCIPstrToRealValue(pos, &(consdata->val[consdata->nvars - 1][consdata->nvarnonz[consdata->nvars - 1]]), &pos);
4303  *success = *success && parsesuccess;
4304  if ( *pos != '\0' )
4305  pos ++; /* remove the "," */
4306 
4307  /* if we got an entry in the upper triangular part, switch the entries for lower triangular */
4308  if ( consdata->col[consdata->nvars - 1][consdata->nvarnonz[consdata->nvars - 1]] >
4309  consdata->row[consdata->nvars - 1][consdata->nvarnonz[consdata->nvars - 1]] )
4310  {
4311  i = consdata->col[consdata->nvars - 1][consdata->nvarnonz[consdata->nvars - 1]];
4312  consdata->col[consdata->nvars - 1][consdata->nvarnonz[consdata->nvars - 1]] =
4313  consdata->row[consdata->nvars - 1][consdata->nvarnonz[consdata->nvars - 1]];
4314  consdata->row[consdata->nvars - 1][consdata->nvarnonz[consdata->nvars - 1]] = i;
4315  }
4316 
4317  consdata->nvarnonz[consdata->nvars - 1]++;
4318 
4319  /* skip whitespace */
4320  while( isspace((unsigned char)*pos) )
4321  pos++;
4322  }
4323 
4324  /* resize the arrays to their final size */
4325  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->col[consdata->nvars - 1], currentsize, consdata->nvarnonz[consdata->nvars - 1]) );
4326  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->row[consdata->nvars - 1], currentsize, consdata->nvarnonz[consdata->nvars - 1]) );
4327  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->val[consdata->nvars - 1], currentsize, consdata->nvarnonz[consdata->nvars - 1]) );
4328 
4329  /* skip whitespace */
4330  while ( isspace((unsigned char)*pos) )
4331  pos++;
4332  }
4333 
4334  /* resize the arrays to their final size */
4335  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->nvarnonz, nvars, consdata->nvars) );
4336  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->col, nvars, consdata->nvars) );
4337  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->row, nvars, consdata->nvars) );
4338  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->val, nvars, consdata->nvars) );
4339  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->vars, nvars, consdata->nvars));
4340 
4341  /* compute sdpnnonz */
4342  for (v = 0; v < consdata->nvars; v++)
4343  consdata->nnonz += consdata->nvarnonz[v];
4344 
4345  /* set maxevsubmat */
4346  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->maxevsubmat, 2) );
4347  consdata->maxevsubmat[0] = -1;
4348  consdata->maxevsubmat[1] = -1;
4349 
4350  /* create the constraint */
4351  SCIP_CALL( SCIPcreateCons(scip, cons, name, conshdlr, consdata, initial, separate, enforce, check, propagate, local, modifiable,
4352  dynamic, removable, stickingatnode) );
4353 
4354  /* compute maximum rhs entry for later use in the DIMACS Error Norm */
4355  SCIP_CALL( setMaxRhsEntry(*cons) );
4356 
4357 #ifdef SCIP_MORE_DEBUG
4358  SCIP_CALL( SCIPprintCons(scip, *cons, NULL) );
4359 #endif
4360 
4361  return SCIP_OKAY;
4362 }
4363 
4365 static
4366 SCIP_DECL_CONSGETVARS(consGetVarsSdp)
4367 {/*lint --e{715}*/
4368  SCIP_CONSDATA* consdata;
4369  int nvars;
4370  int i;
4371 
4372  assert( scip != NULL );
4373  assert( cons != NULL );
4374  assert( vars != NULL );
4375  assert( success != NULL );
4376  assert( varssize >= 0 );
4377 
4378  consdata = SCIPconsGetData(cons);
4379  assert( consdata != NULL );
4380 
4381  nvars = consdata->nvars;
4382 
4383  if ( nvars > varssize )
4384  {
4385  SCIPdebugMsg(scip, "consGetVarsIndicator called for array of size %d, needed size %d.\n", varssize, nvars);
4386  *success = FALSE;
4387  return SCIP_OKAY;
4388  }
4389 
4390  for (i = 0; i < nvars; i++)
4391  vars[i] = consdata->vars[i];
4392 
4393  *success = TRUE;
4394 
4395  return SCIP_OKAY;
4396 }
4397 
4399 static
4400 SCIP_DECL_CONSGETNVARS(consGetNVarsSdp)
4401 {/*lint --e{715}*/
4402  SCIP_CONSDATA* consdata;
4403 
4404  assert( scip != NULL );
4405  assert( cons != NULL );
4406  assert( nvars != NULL );
4407  assert( success != NULL );
4408 
4409  consdata = SCIPconsGetData(cons);
4410  assert( consdata != NULL );
4411 
4412  *nvars = consdata->nvars;
4413  *success = TRUE;
4414 
4415  return SCIP_OKAY;
4416 }
4417 
4419 SCIP_RETCODE SCIPincludeConshdlrSdp(
4420  SCIP* scip
4421  )
4422 {
4423  SCIP_CONSHDLR* conshdlr = NULL;
4424  SCIP_CONSHDLRDATA* conshdlrdata = NULL;
4425 
4426  assert( scip != NULL );
4427 
4428  /* allocate memory for the conshdlrdata */
4429  SCIP_CALL( SCIPallocMemory(scip, &conshdlrdata) );
4430  conshdlrdata->quadconsrank1 = FALSE;
4431  conshdlrdata->upgradquadconss = FALSE;
4432  conshdlrdata->upgradekeepquad = FALSE;
4433  conshdlrdata->quadconsidx = NULL;
4434  conshdlrdata->quadconsvars = NULL;
4435  conshdlrdata->nquadconsidx = 0;
4436  conshdlrdata->X = NULL;
4437  conshdlrdata->nsdpvars = 0;
4438  conshdlrdata->sdpcons = NULL;
4439  conshdlrdata->triedlinearconss = FALSE;
4440  conshdlrdata->sdpconshdlrdata = conshdlrdata; /* set this to itself to simplify access of parameters */
4441 
4442  /* include constraint handler */
4443  SCIP_CALL( SCIPincludeConshdlrBasic(scip, &conshdlr, CONSHDLR_NAME, CONSHDLR_DESC,
4444  CONSHDLR_ENFOPRIORITY, CONSHDLR_CHECKPRIORITY, CONSHDLR_EAGERFREQ, CONSHDLR_NEEDSCONS,
4445  consEnfolpSdp, consEnfopsSdp, consCheckSdp, consLockSdp, conshdlrdata) );
4446 
4447  assert( conshdlr != NULL );
4448 
4449  /* set non-fundamental callbacks via specific setter functions */
4450  SCIP_CALL( SCIPsetConshdlrDelete(scip, conshdlr, consDeleteSdp) );
4451  SCIP_CALL( SCIPsetConshdlrFree(scip, conshdlr, consFreeSdp) );
4452  SCIP_CALL( SCIPsetConshdlrCopy(scip, conshdlr, conshdlrCopySdp, consCopySdp) );
4453  SCIP_CALL( SCIPsetConshdlrInitpre(scip, conshdlr,consInitpreSdp) );
4454  SCIP_CALL( SCIPsetConshdlrExit(scip, conshdlr, consExitSdp) );
4455  SCIP_CALL( SCIPsetConshdlrExitpre(scip, conshdlr, consExitpreSdp) );
4456  SCIP_CALL( SCIPsetConshdlrInitsol(scip, conshdlr, consInitsolSdp) );
4457  SCIP_CALL( SCIPsetConshdlrPresol(scip, conshdlr, consPresolSdp, CONSHDLR_MAXPREROUNDS, CONSHDLR_PRESOLTIMING) );
4458  SCIP_CALL( SCIPsetConshdlrSepa(scip, conshdlr, consSepalpSdp, consSepasolSdp, CONSHDLR_SEPAFREQ,
4459  CONSHDLR_SEPAPRIORITY, CONSHDLR_DELAYSEPA) );
4460  SCIP_CALL( SCIPsetConshdlrEnforelax(scip, conshdlr, consEnforelaxSdp) );
4461  SCIP_CALL( SCIPsetConshdlrTrans(scip, conshdlr, consTransSdp) );
4462  SCIP_CALL( SCIPsetConshdlrPrint(scip, conshdlr, consPrintSdp) );
4463  SCIP_CALL( SCIPsetConshdlrParse(scip, conshdlr, consParseSdp) );
4464  SCIP_CALL( SCIPsetConshdlrGetVars(scip, conshdlr, consGetVarsSdp) );
4465  SCIP_CALL( SCIPsetConshdlrGetNVars(scip, conshdlr, consGetNVarsSdp) );
4466 
4467  /* add parameter */
4468 #ifdef OMP
4469  SCIP_CALL( SCIPaddIntParam(scip, "constraints/SDP/threads", "number of threads used for OpenBLAS",
4470  &(conshdlrdata->nthreads), TRUE, DEFAULT_NTHREADS, 1, INT_MAX, NULL, NULL) );
4471 #endif
4472  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/SDP/diaggezerocuts",
4473  "Should linear cuts enforcing the non-negativity of diagonal entries of SDP-matrices be added?",
4474  &(conshdlrdata->diaggezerocuts), TRUE, DEFAULT_DIAGGEZEROCUTS, NULL, NULL) );
4475 
4476  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/SDP/diagzeroimplcuts",
4477  "Should linear cuts enforcing the implications of diagonal entries of zero in SDP-matrices be added?",
4478  &(conshdlrdata->diagzeroimplcuts), TRUE, DEFAULT_DIAGZEROIMPLCUTS, NULL, NULL) );
4479 
4480  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/SDP/twominorlinconss",
4481  "Should linear cuts corresponding to 2 by 2 minors be added?",
4482  &(conshdlrdata->twominorlinconss), TRUE, DEFAULT_TWOMINORLINCONSS, NULL, NULL) );
4483 
4484  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/SDP/twominorprodconss",
4485  "Should linear cuts corresponding to products of 2 by 2 minors be added?",
4486  &(conshdlrdata->twominorprodconss), TRUE, DEFAULT_TWOMINORPRODCONSS, NULL, NULL) );
4487 
4488  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/SDP/quadconsrank1",
4489  "Should quadratic cons for 2x2 minors be added in the rank-1 case?",
4490  &(conshdlrdata->quadconsrank1), TRUE, DEFAULT_QUADCONSRANK1, NULL, NULL) );
4491 
4492  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/SDP/upgradquadconss",
4493  "Should quadratic constraints be upgraded to a rank 1 SDP?",
4494  &(conshdlrdata->upgradquadconss), TRUE, DEFAULT_UPGRADQUADCONSS, NULL, NULL) );
4495 
4496  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/SDP/upgradekeepquad",
4497  "Should the quadratic constraints be kept in the problem after upgrading and the corresponding SDP constraint be added without the rank 1 constraint?",
4498  &(conshdlrdata->upgradekeepquad), TRUE, DEFAULT_UPGRADEKEEPQUAD, NULL, NULL) );
4499 
4500  SCIP_CALL( SCIPaddIntParam(scip, "constraints/SDP/maxnvarsquadupgd",
4501  "maximal number of quadratic constraints and appearing variables so that the QUADCONSUPGD is performed",
4502  &(conshdlrdata->maxnvarsquadupgd), TRUE, DEFAULT_MAXNVARSQUADUPGD, 0, INT_MAX, NULL, NULL) );
4503 
4504  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/SDP/rank1approxheur",
4505  "Should the heuristic that computes the best rank-1 approximation for a given solution be executed?",
4506  &(conshdlrdata->rank1approxheur), TRUE, DEFAULT_RANK1APPROXHEUR, NULL, NULL) );
4507 
4508  return SCIP_OKAY;
4509 }
4510 
4512 SCIP_RETCODE SCIPincludeConshdlrSdpRank1(
4513  SCIP* scip
4514  )
4515 {
4516  SCIP_CONSHDLR* conshdlr = NULL;
4517  SCIP_CONSHDLR* sdpconshdlr;
4518  SCIP_CONSHDLRDATA* conshdlrdata = NULL;
4519 
4520  assert( scip != NULL );
4521 
4522  /* allocate memory for the conshdlrdata */
4523  SCIP_CALL( SCIPallocMemory(scip, &conshdlrdata) );
4524 
4525  /* only use one parameter */
4526  conshdlrdata->diaggezerocuts = FALSE;
4527  conshdlrdata->diagzeroimplcuts = FALSE;
4528  conshdlrdata->triedlinearconss = FALSE;
4529  conshdlrdata->quadconsrank1 = FALSE;
4530  conshdlrdata->rank1approxheur = FALSE;
4531  conshdlrdata->maxnvarsquadupgd = 0;
4532 
4533  /* parameters are retrieved through the SDP constraint handler */
4534  sdpconshdlr = SCIPfindConshdlr(scip, CONSHDLR_NAME);
4535  if ( sdpconshdlr == NULL )
4536  {
4537  SCIPerrorMessage("Needs constraint handler <%s> to work.\n", CONSHDLR_NAME);
4538  return SCIP_PLUGINNOTFOUND;
4539  }
4540  conshdlrdata->sdpconshdlrdata = SCIPconshdlrGetData(sdpconshdlr);
4541  assert( conshdlrdata->sdpconshdlrdata != NULL );
4542 
4543  conshdlrdata->quadconsidx = NULL;
4544  conshdlrdata->quadconsvars = NULL;
4545  conshdlrdata->nquadconsidx = 0;
4546  conshdlrdata->X = NULL;
4547  conshdlrdata->nsdpvars = 0;
4548  conshdlrdata->sdpcons = NULL;
4549 
4550  /* include constraint handler */
4551  SCIP_CALL( SCIPincludeConshdlrBasic(scip, &conshdlr, CONSHDLRRANK1_NAME, CONSHDLRRANK1_DESC,
4552  CONSHDLR_ENFOPRIORITY, CONSHDLR_CHECKPRIORITY, CONSHDLR_EAGERFREQ, CONSHDLR_NEEDSCONS,
4553  consEnfolpSdp, consEnfopsSdp, consCheckSdp, consLockSdp, conshdlrdata) );
4554 
4555  assert( conshdlr != NULL );
4556 
4557  /* set non-fundamental callbacks via specific setter functions */
4558  SCIP_CALL( SCIPsetConshdlrDelete(scip, conshdlr, consDeleteSdp) );
4559  SCIP_CALL( SCIPsetConshdlrFree(scip, conshdlr, consFreeSdp) );
4560  SCIP_CALL( SCIPsetConshdlrCopy(scip, conshdlr, conshdlrCopySdpRank1, consCopySdp) );
4561  SCIP_CALL( SCIPsetConshdlrInitpre(scip, conshdlr,consInitpreSdp) );
4562  SCIP_CALL( SCIPsetConshdlrExit(scip, conshdlr, consExitSdp) );
4563  SCIP_CALL( SCIPsetConshdlrExitpre(scip, conshdlr, consExitpreSdp) );
4564  SCIP_CALL( SCIPsetConshdlrInitsol(scip, conshdlr, consInitsolSdp) );
4565  SCIP_CALL( SCIPsetConshdlrPresol(scip, conshdlr, consPresolSdp, CONSHDLR_MAXPREROUNDS, CONSHDLR_PRESOLTIMING) );
4566  SCIP_CALL( SCIPsetConshdlrSepa(scip, conshdlr, consSepalpSdp, consSepasolSdp, CONSHDLR_SEPAFREQ,
4567  CONSHDLR_SEPAPRIORITY, CONSHDLR_DELAYSEPA) );
4568  SCIP_CALL( SCIPsetConshdlrEnforelax(scip, conshdlr, consEnforelaxSdp) );
4569  SCIP_CALL( SCIPsetConshdlrTrans(scip, conshdlr, consTransSdp) );
4570  SCIP_CALL( SCIPsetConshdlrPrint(scip, conshdlr, consPrintSdp) );
4571  SCIP_CALL( SCIPsetConshdlrParse(scip, conshdlr, consParseSdp) );
4572  SCIP_CALL( SCIPsetConshdlrGetVars(scip, conshdlr, consGetVarsSdp) );
4573  SCIP_CALL( SCIPsetConshdlrGetNVars(scip, conshdlr, consGetNVarsSdp) );
4574 
4575  /* include upgrading function for quadratic constraints */
4576  SCIP_CALL( SCIPincludeQuadconsUpgrade(scip, consQuadConsUpgdSdp, 0, TRUE, CONSHDLRRANK1_NAME) );
4577 
4578  return SCIP_OKAY;
4579 }
4580 
4584 int SCIPconsSdpCompLowerTriangPos(
4585  int i,
4586  int j
4587  )
4588 {
4589  assert( j >= 0 );
4590  assert( i >= j );
4591 
4592  return i*(i+1)/2 + j;
4593 }
4594 
4602 SCIP_RETCODE SCIPconsSdpGetData(
4603  SCIP* scip,
4604  SCIP_CONS* cons,
4605  int* nvars,
4606  int* nnonz,
4607  int* blocksize,
4608  int* arraylength,
4609  int* nvarnonz,
4611  int** col,
4612  int** row,
4613  SCIP_Real** val,
4614  SCIP_VAR** vars,
4615  int* constnnonz,
4617  int* constcol,
4618  int* constrow,
4619  SCIP_Real* constval,
4620  SCIP_Bool* rankone,
4621  int** maxevsubmat,
4622  SCIP_Bool* addedquadcons
4623  )
4624 {
4625  SCIP_CONSDATA* consdata;
4626  int i;
4627 
4628  assert( scip != NULL );
4629  assert( cons != NULL );
4630  assert( nvars != NULL );
4631  assert( nnonz != NULL );
4632  assert( blocksize != NULL );
4633  assert( arraylength != NULL );
4634  assert( nvarnonz != NULL );
4635  assert( col != NULL );
4636  assert( row != NULL );
4637  assert( val != NULL );
4638  assert( vars != NULL );
4639  assert( constnnonz != NULL );
4640 
4641  consdata = SCIPconsGetData(cons);
4642 
4643  assert( consdata->constnnonz == 0 || ( constcol != NULL && constrow != NULL && constval != NULL ) );
4644 
4645  *nvars = consdata->nvars;
4646  *nnonz = consdata->nnonz;
4647  *blocksize = consdata->blocksize;
4648 
4649  for (i = 0; i < consdata->nvars; i++)
4650  vars[i] = consdata->vars[i];
4651 
4652  /* check that the sdp-arrays are long enough to store the information */
4653  if ( *arraylength < consdata->nvars )
4654  {
4655  SCIPdebugMsg(scip, "nvarnonz, col, row and val arrays were not long enough to store the information for cons %s, they need to be at least"
4656  "size %d, given was only length %d! \n", SCIPconsGetName(cons), consdata->nvars, *arraylength);
4657  *arraylength = consdata->nvars;
4658  }
4659  else
4660  {
4661  for (i = 0; i < consdata->nvars; i++)
4662  {
4663  nvarnonz[i] = consdata->nvarnonz[i];
4664  /* set the pointers for each variable */
4665  col[i] = consdata->col[i];
4666  row[i] = consdata->row[i];
4667  val[i] = consdata->val[i];
4668  }
4669  }
4670 
4671  /* set the constant pointers (if a constant part exists) */
4672  if ( consdata->constnnonz > 0 )
4673  {
4674  if ( consdata->constnnonz > *constnnonz )
4675  {
4676  SCIPdebugMsg(scip, "The constant nonzeros arrays were not long enough to store the information for cons %s, they need to be at least"
4677  "size %d, given was only length %d! \n", SCIPconsGetName(cons), consdata->constnnonz, *constnnonz);
4678  }
4679  else
4680  {
4681  for (i = 0; i < consdata->constnnonz; i++)
4682  {
4683  constcol[i] = consdata->constcol[i];
4684  constrow[i] = consdata->constrow[i];
4685  constval[i] = consdata->constval[i];
4686  }
4687  }
4688  }
4689 
4690  *constnnonz = consdata->constnnonz;
4691 
4692  /* set the information about rankone, the current submatrix with largest minimal eigenvalue ([-1,-1] if not yet
4693  computed), and the quadratic 2x2-minor constraints if desired */
4694  if ( rankone != NULL && maxevsubmat != NULL )
4695  {
4696  *rankone = consdata->rankone;
4697  *maxevsubmat[0] = consdata->maxevsubmat[0];
4698  *maxevsubmat[1] = consdata->maxevsubmat[1];
4699  *addedquadcons = consdata->addedquadcons;
4700  }
4701 
4702  return SCIP_OKAY;
4703 }
4704 
4709 SCIP_RETCODE SCIPconsSdpGetNNonz(
4710  SCIP* scip,
4711  SCIP_CONS* cons,
4712  int* nnonz,
4713  int* constnnonz
4714  )
4715 {
4716  SCIP_CONSDATA* consdata;
4717 
4718  assert( scip != NULL );
4719  assert( cons != NULL );
4720 
4721  consdata = SCIPconsGetData(cons);
4722  assert( consdata != NULL );
4723 
4724  if ( nnonz != NULL )
4725  *nnonz = consdata->nnonz;
4726 
4727  if ( constnnonz != NULL )
4728  *constnnonz = consdata->constnnonz;
4729 
4730  return SCIP_OKAY;
4731 }
4732 
4734 int SCIPconsSdpGetBlocksize(
4735  SCIP* scip,
4736  SCIP_CONS* cons
4737  )
4738 {
4739  SCIP_CONSDATA* consdata;
4740 
4741  assert( scip != NULL );
4742  assert( cons != NULL );
4743 
4744  consdata = SCIPconsGetData(cons);
4745  assert( consdata != NULL );
4746 
4747  return consdata->blocksize;
4748 }
4749 
4751 SCIP_RETCODE SCIPconsSdpGetFullAj(
4752  SCIP* scip,
4753  SCIP_CONS* cons,
4754  int j,
4755  SCIP_Real* Aj
4756  )
4757 {
4758  SCIP_CONSDATA* consdata;
4759  int blocksize;
4760  int i;
4761 
4762  assert( scip != NULL );
4763  assert( cons != NULL );
4764  assert( j >= 0 );
4765  assert( Aj != NULL );
4766 
4767  consdata = SCIPconsGetData(cons);
4768  assert( consdata != NULL );
4769  blocksize = consdata->blocksize;
4770 
4771  assert( j < consdata->nvars );
4772 
4773  for (i = 0; i < blocksize * blocksize; i++)
4774  Aj[i] = 0;
4775 
4776  for (i = 0; i < consdata->nvarnonz[j]; i++)
4777  {
4778  Aj[consdata->col[j][i] * blocksize + consdata->row[j][i]] = consdata->val[j][i]; /*lint !e679*/
4779  Aj[consdata->row[j][i] * blocksize + consdata->col[j][i]] = consdata->val[j][i]; /*lint !e679*/
4780  }
4781 
4782  return SCIP_OKAY;
4783 }
4784 
4786 SCIP_RETCODE SCIPconsSdpGetFullConstMatrix(
4787  SCIP* scip,
4788  SCIP_CONS* cons,
4789  SCIP_Real* mat
4790  )
4791 {
4792  SCIP_CONSDATA* consdata;
4793  int blocksize;
4794  int i;
4795  int j;
4796 
4797  assert( scip != NULL );
4798  assert( cons != NULL );
4799  assert( mat != NULL );
4800 
4801  consdata = SCIPconsGetData(cons);
4802  blocksize = consdata->blocksize;
4803 
4804  for (i = 0; i < blocksize; i++)
4805  {
4806  for (j = 0; j < blocksize; j++)
4807  mat[i * blocksize + j] = 0.0; /*lint !e679*/
4808  }
4809 
4810  for (i = 0; i < consdata->constnnonz; i++)
4811  {
4812  mat[consdata->constcol[i] * blocksize + consdata->constrow[i]] = consdata->constval[i]; /*lint !e679*/
4813  mat[consdata->constrow[i] * blocksize + consdata->constcol[i]] = consdata->constval[i]; /*lint !e679*/
4814  }
4815 
4816  return SCIP_OKAY;
4817 }
4818 
4820 SCIP_RETCODE SCIPconsSdpGetLowerTriangConstMatrix(
4821  SCIP* scip,
4822  SCIP_CONS* cons,
4823  SCIP_Real* mat
4824  )
4825 {
4826  SCIP_CONSDATA* consdata;
4827  int blocksize;
4828  int i;
4829 
4830  assert( scip != NULL );
4831  assert( cons != NULL );
4832  assert( mat != NULL );
4833 
4834  consdata = SCIPconsGetData(cons);
4835  assert( consdata != NULL );
4836 
4837  blocksize = consdata->blocksize;
4838 
4839  /* initialize the matrix with 0 */
4840  for (i = 0; i < (blocksize * (blocksize + 1)) / 2; i++)
4841  mat[i] = 0.0;
4842 
4843  for (i = 0; i < consdata->constnnonz; i++)
4844  mat[SCIPconsSdpCompLowerTriangPos(consdata->constrow[i], consdata->constcol[i])] = consdata->constval[i];
4845 
4846  return SCIP_OKAY;
4847 }
4848 
4859 SCIP_RETCODE SCIPconsSdpGuessInitialPoint(
4860  SCIP* scip,
4861  SCIP_CONS* cons,
4862  SCIP_Real* lambdastar
4863  )
4864 {
4865  SCIP_CONSDATA* consdata;
4866  SCIP_Real sparsity;
4867  SCIP_Real maxinfnorm;
4868  SCIP_Real maxconst;
4869  SCIP_Real mininfnorm;
4870  SCIP_Real maxobj;
4871  SCIP_Real maxbound;
4872  SCIP_Real primalguess;
4873  SCIP_Real dualguess;
4874  SCIP_Real compval;
4875  int blocksize;
4876  int i;
4877  int v;
4878 
4879  assert( scip != NULL );
4880  assert( cons != NULL );
4881  assert( lambdastar != NULL );
4882  assert( strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(cons)), CONSHDLR_NAME) == 0 || strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(cons)), CONSHDLRRANK1_NAME) == 0 );
4883 
4884  consdata = SCIPconsGetData(cons);
4885  assert( consdata != NULL );
4886 
4887  /* If there are no nonzeros, we cannot use the usual formula, since it divides through the number of nonzeros. In this case,
4888  * however, we will not solve an SDP anyways but at most an LP (more likely the problem will be solved in local presolving,
4889  * if all variables are fixed and not only those in the SDP-part), so we just take the default value of SDPA.
4890  */
4891  if ( consdata->nnonz == 0 )
4892  {
4893  *lambdastar = 100.0;
4894  return SCIP_OKAY;
4895  }
4896 
4897  blocksize = consdata->blocksize;
4898 
4899  sparsity = consdata->nnonz / (0.5 * blocksize * (blocksize + 1));
4900 
4901  /* compute the maximum entry of the A_i */
4902  maxinfnorm = 0.0;
4903  mininfnorm = SCIPinfinity(scip);
4904  for (v = 0; v < consdata->nvars; v++)
4905  {
4906  for (i = 0; i < consdata->nvarnonz[v]; i++)
4907  {
4908  if ( SCIPisGT(scip, REALABS(consdata->val[v][i]), maxinfnorm ) )
4909  maxinfnorm = REALABS(consdata->val[v][i]);
4910  if ( SCIPisLT(scip, REALABS(consdata->val[v][i]), mininfnorm) )
4911  mininfnorm = REALABS(consdata->val[v][i]);
4912  }
4913  }
4914  maxconst = 0.0;
4915  for (i = 0; i < consdata->constnnonz; i++)
4916  {
4917  if ( SCIPisGT(scip, REALABS(consdata->constval[i]), maxconst ) )
4918  maxconst = REALABS(consdata->constval[i]);
4919  }
4920 
4921  assert( SCIPisGT(scip, mininfnorm, 0.0) );
4922 
4923  /* compute maximum b_i and bound */
4924  maxobj = 0.0;
4925  maxbound = 0.0;
4926  for (v = 0; v < consdata->nvars; v++)
4927  {
4928  if ( SCIPisGT(scip, REALABS(SCIPvarGetObj(consdata->vars[v])), maxobj) )
4929  maxobj = REALABS(SCIPvarGetObj(consdata->vars[v]));
4930  compval = SCIPisInfinity(scip, REALABS(SCIPvarGetUbGlobal(consdata->vars[v]))) ? 1e+6 : REALABS(SCIPvarGetUbGlobal(consdata->vars[v]));
4931  if ( SCIPisGT(scip, compval, maxbound) )
4932  maxbound = compval;
4933  compval = SCIPisInfinity(scip, REALABS(SCIPvarGetLbGlobal(consdata->vars[v]))) ? 1e+6 : REALABS(SCIPvarGetUbGlobal(consdata->vars[v]));
4934  if ( SCIPisGT(scip, compval, maxbound) )
4935  maxbound = compval;
4936  }
4937 
4938  /* if all variables were unbounded, we set the value to 10^6 */
4939  if ( SCIPisEQ(scip, maxbound, 0.0) )
4940  maxbound = 1E+6;
4941 
4942  /* compute primal and dual guess */
4943  primalguess = maxobj / (sparsity * mininfnorm);
4944  dualguess = sparsity * maxinfnorm * maxbound + maxconst;
4945 
4946  if ( SCIPisGT(scip, primalguess, dualguess) )
4947  *lambdastar = primalguess;
4948  else
4949  *lambdastar = dualguess;
4950 
4951  return SCIP_OKAY;
4952 }
4953 
4955 SCIP_Real SCIPconsSdpGetMaxConstEntry(
4956  SCIP* scip,
4957  SCIP_CONS* cons
4958  )
4959 {
4960  SCIP_CONSDATA* consdata;
4961 
4962  assert( scip != NULL );
4963  assert( cons != NULL );
4964 
4965  consdata = SCIPconsGetData(cons);
4966 
4967  return consdata->maxrhsentry;
4968 }
4969 
4971 SCIP_Real SCIPconsSdpGetMaxSdpCoef(
4972  SCIP* scip,
4973  SCIP_CONS* cons
4974  )
4975 {
4976  SCIP_CONSDATA* consdata;
4977  SCIP_Real maxcoef;
4978  int v;
4979  int i;
4980 
4981  assert( scip != NULL );
4982  assert( cons != NULL );
4983 
4984  consdata = SCIPconsGetData(cons);
4985 
4986  maxcoef = 0.0;
4987 
4988  for (v = 0; v < consdata->nvars; v++)
4989  {
4990  for (i = 0; i < consdata->nvarnonz[v]; i++)
4991  {
4992  if ( SCIPisGT(scip, REALABS(consdata->val[v][i]), maxcoef) )
4993  maxcoef = REALABS(consdata->val[v][i]);
4994  }
4995  }
4996 
4997  return maxcoef;
4998 }
4999 
5005 int SCIPconsSdpComputeUbSparseSdpMatrixLength(
5006  SCIP_CONS* cons
5007  )
5008 { /*lint --e{715}*/
5009  SCIP_CONSDATA* consdata;
5010  int v;
5011  int ub;
5012  int denselength;
5013 
5014  assert( cons != NULL );
5015 
5016  consdata = SCIPconsGetData(cons);
5017  assert( consdata != NULL );
5018 
5019  ub = consdata->constnnonz;
5020 
5021  for (v = 0; v < consdata->nvars; v++)
5022  ub += consdata->nvarnonz[v];
5023 
5024  denselength = consdata->blocksize * (consdata->blocksize + 1) / 2;
5025 
5026  return (ub <= denselength ? ub : denselength);
5027 }
5028 
5033 SCIP_RETCODE SCIPconsSdpComputeSparseSdpMatrix(
5034  SCIP* scip,
5035  SCIP_CONS* cons,
5036  SCIP_SOL* sol,
5037  int* length,
5039  int* row,
5040  int* col,
5041  SCIP_Real* val
5042  )
5043 {
5044  SCIP_CONSDATA* consdata;
5045  int i;
5046  int v;
5047  int nnonz;
5048 
5049  assert( scip != NULL );
5050  assert( cons != NULL );
5051  assert( sol != NULL );
5052  assert( length != NULL );
5053  assert( row != NULL );
5054  assert( col != NULL );
5055  assert( val != NULL );
5056 
5057  consdata = SCIPconsGetData(cons);
5058  assert( consdata != NULL );
5059 
5060  /* initialize nnonz/row/col/val with constant arrays */
5061  nnonz = consdata->constnnonz;
5062  if ( *length < nnonz )
5063  {
5064  *length = -1;
5065  SCIPerrorMessage("Arrays not long enough in SCIPconsSdpComputeSparseSdpMatrix, length %d given, need at least %d (probably more)\n",
5066  *length, nnonz);
5067  return SCIP_ERROR;
5068  }
5069 
5070  for (i = 0; i < consdata->constnnonz; i++)
5071  {
5072  row[i] = consdata->constrow[i];
5073  col[i] = consdata->constcol[i];
5074  val[i] = -1 * consdata->constval[i];
5075  }
5076 
5077  /* add all variable arrays multiplied by corresponding solution value */
5078  for (v = 0; v < consdata->nvars; v++)
5079  {
5080  SCIP_CALL( SCIPsdpVarfixerMergeArrays(SCIPblkmem(scip), SCIPepsilon(scip), consdata->row[v], consdata->col[v], consdata->val[v], consdata->nvarnonz[v],
5081  FALSE, SCIPgetSolVal(scip, sol, consdata->vars[v]), row, col, val, &nnonz, *length) );
5082  if ( nnonz > *length )
5083  {
5084  *length = -1;
5085  SCIPerrorMessage("Arrays not long enough in SCIPconsSdpComputeSparseSdpMatrix, length %d given, need at least %d (probably more)\n",
5086  *length, nnonz);
5087  return SCIP_ERROR;
5088  }
5089  }
5090 
5091  /* update length pointer */
5092  *length = nnonz;
5093 
5094  return SCIP_OKAY;
5095 }
5096 
5098 SCIP_Bool SCIPconsSdpShouldBeRankOne(
5099  SCIP_CONS* cons
5100  )
5101 {
5102  SCIP_CONSDATA* consdata;
5103 
5104  assert( cons != NULL );
5105 
5106  consdata = SCIPconsGetData(cons);
5107  assert( consdata != NULL );
5108 
5109  return consdata->rankone;
5110 }
5111 
5113 SCIP_RETCODE SCIPconsSdpGetMaxEVSubmat(
5114  SCIP_CONS* cons,
5115  int** maxevsubmat
5117  )
5118 {
5119  SCIP_CONSDATA* consdata;
5120 
5121  assert( cons != NULL );
5122 
5123  consdata = SCIPconsGetData(cons);
5124  assert( consdata != NULL );
5125  assert( maxevsubmat != NULL );
5126 
5127  *maxevsubmat[0] = consdata->maxevsubmat[0];
5128  *maxevsubmat[1] = consdata->maxevsubmat[1];
5129 
5130  return SCIP_OKAY;
5131 }
5132 
5134 SCIP_Bool SCIPconsSdpAddedQuadCons(
5135  SCIP_CONS* cons
5136  )
5137 {
5138  SCIP_CONSDATA* consdata;
5139 
5140  assert( cons != NULL );
5141 
5142  consdata = SCIPconsGetData(cons);
5143  assert( consdata != NULL );
5144 
5145  return consdata->addedquadcons;
5146 }
5147 
5149 SCIP_RETCODE SCIPcreateConsSdp(
5150  SCIP* scip,
5151  SCIP_CONS** cons,
5152  const char* name,
5153  int nvars,
5154  int nnonz,
5155  int blocksize,
5156  int* nvarnonz,
5157  int** col,
5158  int** row,
5159  SCIP_Real** val,
5160  SCIP_VAR** vars,
5161  int constnnonz,
5162  int* constcol,
5163  int* constrow,
5164  SCIP_Real* constval
5165  )
5166 {
5167  SCIP_CONSHDLR* conshdlr;
5168  SCIP_CONSDATA* consdata = NULL;
5169  int i;
5170  int j;
5171 
5172  assert( scip != NULL );
5173  assert( cons != NULL );
5174  assert( name != NULL );
5175  assert( nvars >= 0 );
5176  assert( nnonz >= 0 );
5177  assert( blocksize >= 0 );
5178  assert( constnnonz >= 0 );
5179  assert( nvars == 0 || vars != NULL );
5180  assert( nnonz == 0 || (nvarnonz != NULL && col != NULL && row != NULL && val != NULL ));
5181  assert( constnnonz == 0 || (constcol != NULL && constrow != NULL && constval != NULL ));
5182 
5183  conshdlr = SCIPfindConshdlr(scip, "SDP");
5184  if ( conshdlr == NULL )
5185  {
5186  SCIPerrorMessage("SDP constraint handler not found\n");
5187  return SCIP_PLUGINNOTFOUND;
5188  }
5189 
5190  /* create constraint data */
5191  SCIP_CALL( SCIPallocBlockMemory(scip, &consdata) );
5192  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->nvarnonz, nvars) );
5193  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->col, nvars) );
5194  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->row, nvars) );
5195  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->val, nvars) );
5196  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->constcol, constnnonz) );
5197  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->constrow, constnnonz) );
5198  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->constval, constnnonz) );
5199  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->vars, nvars) );
5200 
5201  for (i = 0; i < nvars; i++)
5202  {
5203  assert( nvarnonz[i] >= 0 );
5204 
5205  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->col[i], nvarnonz[i]));
5206  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->row[i], nvarnonz[i]));
5207  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->val[i], nvarnonz[i]));
5208  }
5209 
5210  consdata->nvars = nvars;
5211  consdata->nnonz = nnonz;
5212  consdata->constnnonz = constnnonz;
5213  consdata->blocksize = blocksize;
5214  consdata->locks = NULL;
5215 
5216  for (i = 0; i < nvars; i++)
5217  {
5218  consdata->nvarnonz[i] = nvarnonz[i];
5219 
5220  for (j = 0; j < nvarnonz[i]; j++)
5221  {
5222  assert( col[i][j] >= 0 );
5223  assert( col[i][j] < blocksize );
5224  assert( row[i][j] >= 0 );
5225  assert( row[i][j] < blocksize );
5226 
5227  consdata->col[i][j] = col[i][j];
5228  consdata->row[i][j] = row[i][j];
5229  consdata->val[i][j] = val[i][j];
5230  }
5231  }
5232 
5233  for (i = 0; i < constnnonz; i++)
5234  {
5235  consdata->constcol[i] = constcol[i];
5236  consdata->constrow[i] = constrow[i];
5237  consdata->constval[i] = constval[i];
5238  }
5239 
5240  for (i = 0; i < nvars; i++)
5241  {
5242  consdata->vars[i] = vars[i];
5243  SCIP_CALL( SCIPcaptureVar(scip, consdata->vars[i]) );
5244  }
5245  SCIPdebugMsg(scip, "creating cons %s.\n", name);
5246 
5247  consdata->rankone = FALSE;
5248 
5249  /* allocate memory for rank one constraint */
5250  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->maxevsubmat, 2) );
5251  consdata->maxevsubmat[0] = -1;
5252  consdata->maxevsubmat[1] = -1;
5253 
5254  /* quadratic 2x2-minor constraints added? */
5255  consdata->addedquadcons = FALSE;
5256 
5257  /* create constraint */
5258  SCIP_CALL( SCIPcreateCons(scip, cons, name, conshdlr, consdata, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE) );
5259 
5260  /* compute maximum rhs entry for later use in the DIMACS Error Norm */
5261  SCIP_CALL( setMaxRhsEntry(*cons) );
5262 
5263  return SCIP_OKAY;
5264 }
5265 
5266 
5268 SCIP_RETCODE SCIPcreateConsSdpRank1(
5269  SCIP* scip,
5270  SCIP_CONS** cons,
5271  const char* name,
5272  int nvars,
5273  int nnonz,
5274  int blocksize,
5275  int* nvarnonz,
5276  int** col,
5277  int** row,
5278  SCIP_Real** val,
5279  SCIP_VAR** vars,
5280  int constnnonz,
5281  int* constcol,
5282  int* constrow,
5283  SCIP_Real* constval
5284  )
5285 {
5286  SCIP_CONSHDLR* conshdlr;
5287  SCIP_CONSDATA* consdata = NULL;
5288  int i;
5289  int j;
5290 
5291  assert( scip != NULL );
5292  assert( cons != NULL );
5293  assert( name != NULL );
5294  assert( nvars >= 0 );
5295  assert( nnonz >= 0 );
5296  assert( blocksize >= 0 );
5297  assert( constnnonz >= 0 );
5298  assert( nvars == 0 || vars != NULL );
5299  assert( nnonz == 0 || (nvarnonz != NULL && col != NULL && row != NULL && val != NULL ));
5300  assert( constnnonz == 0 || (constcol != NULL && constrow != NULL && constval != NULL ));
5301 
5302  conshdlr = SCIPfindConshdlr(scip, CONSHDLRRANK1_NAME);
5303  if ( conshdlr == NULL )
5304  {
5305  SCIPerrorMessage("Rank 1 SDP constraint handler not found\n");
5306  return SCIP_PLUGINNOTFOUND;
5307  }
5308 
5309  /* create constraint data */
5310  SCIP_CALL( SCIPallocBlockMemory(scip, &consdata) );
5311  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->nvarnonz, nvars) );
5312  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->col, nvars) );
5313  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->row, nvars) );
5314  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->val, nvars) );
5315  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->constcol, constnnonz) );
5316  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->constrow, constnnonz) );
5317  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->constval, constnnonz) );
5318  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->vars, nvars) );
5319 
5320  for (i = 0; i < nvars; i++)
5321  {
5322  assert( nvarnonz[i] >= 0 );
5323 
5324  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->col[i], nvarnonz[i]));
5325  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->row[i], nvarnonz[i]));
5326  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->val[i], nvarnonz[i]));
5327  }
5328 
5329  consdata->nvars = nvars;
5330  consdata->nnonz = nnonz;
5331  consdata->constnnonz = constnnonz;
5332  consdata->blocksize = blocksize;
5333  consdata->locks = NULL;
5334 
5335  for (i = 0; i < nvars; i++)
5336  {
5337  consdata->nvarnonz[i] = nvarnonz[i];
5338 
5339  for (j = 0; j < nvarnonz[i]; j++)
5340  {
5341  assert( col[i][j] >= 0 );
5342  assert( col[i][j] < blocksize );
5343  assert( row[i][j] >= 0 );
5344  assert( row[i][j] < blocksize );
5345 
5346  consdata->col[i][j] = col[i][j];
5347  consdata->row[i][j] = row[i][j];
5348  consdata->val[i][j] = val[i][j];
5349  }
5350  }
5351 
5352  for (i = 0; i < constnnonz; i++)
5353  {
5354  consdata->constcol[i] = constcol[i];
5355  consdata->constrow[i] = constrow[i];
5356  consdata->constval[i] = constval[i];
5357  }
5358 
5359  for (i = 0; i < nvars; i++)
5360  {
5361  consdata->vars[i] = vars[i];
5362  SCIP_CALL( SCIPcaptureVar(scip, consdata->vars[i]) );
5363  }
5364  SCIPdebugMsg(scip, "creating cons %s (rank 1).\n", name);
5365  consdata->rankone = TRUE;
5366 
5367  /* allocate memory for rank one constraint */
5368  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->maxevsubmat, 2) );
5369  consdata->maxevsubmat[0] = -1;
5370  consdata->maxevsubmat[1] = -1;
5371 
5372  /* quadratic 2x2-minor constraints added? */
5373  consdata->addedquadcons = FALSE;
5374 
5375  /* create constraint */
5376  SCIP_CALL( SCIPcreateCons(scip, cons, name, conshdlr, consdata, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE) );
5377 
5378  /* compute maximum rhs entry for later use in the DIMACS Error Norm */
5379  SCIP_CALL( setMaxRhsEntry(*cons) );
5380 
5381  return SCIP_OKAY;
5382 }
SCIPconsSdpGetMaxConstEntry
SCIP_Real SCIPconsSdpGetMaxConstEntry(SCIP *scip, SCIP_CONS *cons)
Definition: cons_sdp.c:4956
SCIP_DECL_CONSTRANS
static SCIP_DECL_CONSTRANS(consTransSdp)
Definition: cons_sdp.c:3200
addTwoMinorProdConstraints
static SCIP_RETCODE addTwoMinorProdConstraints(SCIP *scip, SCIP_CONS **conss, int nconss, int *naddconss)
Definition: cons_sdp.c:1275
SCIPconsSdpComputeUbSparseSdpMatrixLength
int SCIPconsSdpComputeUbSparseSdpMatrixLength(SCIP_CONS *cons)
Definition: cons_sdp.c:5006
separateSol
static SCIP_RETCODE separateSol(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_CONS *cons, SCIP_SOL *sol, SCIP_RESULT *result)
Definition: cons_sdp.c:603
DEFAULT_QUADCONSRANK1
#define DEFAULT_QUADCONSRANK1
Definition: cons_sdp.c:102
SCIPconsSdpGetFullAj
SCIP_RETCODE SCIPconsSdpGetFullAj(SCIP *scip, SCIP_CONS *cons, int j, SCIP_Real *Aj)
Definition: cons_sdp.c:4752
SCIPcreateConsSdpRank1
SCIP_RETCODE SCIPcreateConsSdpRank1(SCIP *scip, SCIP_CONS **cons, const char *name, int nvars, int nnonz, int blocksize, int *nvarnonz, int **col, int **row, SCIP_Real **val, SCIP_VAR **vars, int constnnonz, int *constcol, int *constrow, SCIP_Real *constval)
Definition: cons_sdp.c:5269
SCIPsdpVarfixerSortRowCol
void SCIPsdpVarfixerSortRowCol(int *row, int *col, SCIP_Real *val, int length)
Definition: SdpVarfixer.c:50
convertRowToColFormatFullMatrix
static SCIP_RETCODE convertRowToColFormatFullMatrix(int rows, int cols, SCIP_Real *rowmatrix, SCIP_Real *colmatrix)
Definition: cons_sdp.c:165
CONSHDLR_PRESOLTIMING
#define CONSHDLR_PRESOLTIMING
Definition: cons_sdp.c:95
CONSHDLR_NAME
#define CONSHDLR_NAME
Definition: cons_sdp.c:78
SCIPincludeConshdlrSdp
SCIP_RETCODE SCIPincludeConshdlrSdp(SCIP *scip)
Definition: cons_sdp.c:4420
multiplyConstraintMatrix
static SCIP_RETCODE multiplyConstraintMatrix(SCIP_CONS *cons, int j, SCIP_Real *v, SCIP_Real *vAv)
Definition: cons_sdp.c:398
PARSE_SIZEFACTOR
#define PARSE_SIZEFACTOR
Definition: cons_sdp.c:97
SCIP_DECL_CONSENFORELAX
static SCIP_DECL_CONSENFORELAX(consEnforelaxSdp)
Definition: cons_sdp.c:3791
enforceConstraint
static SCIP_RETCODE enforceConstraint(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_CONS **conss, int nconss, SCIP_SOL *sol, SCIP_RESULT *result)
Definition: cons_sdp.c:2191
CONSHDLR_NEEDSCONS
#define CONSHDLR_NEEDSCONS
Definition: cons_sdp.c:93
DEFAULT_TWOMINORPRODCONSS
#define DEFAULT_TWOMINORPRODCONSS
Definition: cons_sdp.c:101
SCIP_DECL_CONSDELETE
static SCIP_DECL_CONSDELETE(consDeleteSdp)
Definition: cons_sdp.c:3837
SCIPconsSdpAddedQuadCons
SCIP_Bool SCIPconsSdpAddedQuadCons(SCIP_CONS *cons)
Definition: cons_sdp.c:5135
CONSHDLR_SEPAFREQ
#define CONSHDLR_SEPAFREQ
Definition: cons_sdp.c:87
SCIP_DECL_CONSEXIT
static SCIP_DECL_CONSEXIT(consExitSdp)
Definition: cons_sdp.c:2850
CONSHDLR_DELAYSEPA
#define CONSHDLR_DELAYSEPA
Definition: cons_sdp.c:92
computeSdpMatrix
static SCIP_RETCODE computeSdpMatrix(SCIP *scip, SCIP_CONS *cons, SCIP_SOL *y, SCIP_Real *matrix)
Definition: cons_sdp.c:256
checkVarsLocks
static SCIP_RETCODE checkVarsLocks(SCIP *scip, SCIP_CONS *cons)
Definition: cons_sdp.c:1682
CONSHDLR_EAGERFREQ
#define CONSHDLR_EAGERFREQ
Definition: cons_sdp.c:88
DEFAULT_MAXNVARSQUADUPGD
#define DEFAULT_MAXNVARSQUADUPGD
Definition: cons_sdp.c:105
DEFAULT_DIAGGEZEROCUTS
#define DEFAULT_DIAGGEZEROCUTS
Definition: cons_sdp.c:98
SCIPconsSdpComputeSparseSdpMatrix
SCIP_RETCODE SCIPconsSdpComputeSparseSdpMatrix(SCIP *scip, SCIP_CONS *cons, SCIP_SOL *sol, int *length, int *row, int *col, SCIP_Real *val)
Definition: cons_sdp.c:5034
SCIPconsSdpGetLowerTriangConstMatrix
SCIP_RETCODE SCIPconsSdpGetLowerTriangConstMatrix(SCIP *scip, SCIP_CONS *cons, SCIP_Real *mat)
Definition: cons_sdp.c:4821
SCIP_DECL_CONSSEPALP
static SCIP_DECL_CONSSEPALP(consSepalpSdp)
Definition: cons_sdp.c:3820
setMaxRhsEntry
static SCIP_RETCODE setMaxRhsEntry(SCIP_CONS *cons)
Definition: cons_sdp.c:438
move_1x1_blocks_to_lp
static SCIP_RETCODE move_1x1_blocks_to_lp(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_CONS **conss, int nconss, int *naddconss, int *ndelconss, int *nchgbds, SCIP_RESULT *result)
Definition: cons_sdp.c:1415
SCIPlapackMatrixMatrixMult
SCIP_RETCODE SCIPlapackMatrixMatrixMult(int nrowsA, int ncolsA, SCIP_Real *matrixA, SCIP_Bool transposeA, int nrowsB, int ncolsB, SCIP_Real *matrixB, SCIP_Bool transposeB, SCIP_Real *result)
Definition: lapack_interface.c:393
DEFAULT_RANK1APPROXHEUR
#define DEFAULT_RANK1APPROXHEUR
Definition: cons_sdp.c:106
SCIP_DECL_CONSCHECK
static SCIP_DECL_CONSCHECK(consCheckSdp)
Definition: cons_sdp.c:3307
CONSHDLR_DESC
#define CONSHDLR_DESC
Definition: cons_sdp.c:79
cons_sdp.h
Constraint handler for SDP-constraints.
expandSymMatrix
static SCIP_RETCODE expandSymMatrix(int size, SCIP_Real *symMat, SCIP_Real *fullMat)
Definition: cons_sdp.c:223
SdpVarmapper.h
maps SCIP variables to SDP indices (the SCIP variables are given SDP indices in the order in which th...
SCIPconsSdpGetFullConstMatrix
SCIP_RETCODE SCIPconsSdpGetFullConstMatrix(SCIP *scip, SCIP_CONS *cons, SCIP_Real *mat)
Definition: cons_sdp.c:4787
SCIPcreateConsSdp
SCIP_RETCODE SCIPcreateConsSdp(SCIP *scip, SCIP_CONS **cons, const char *name, int nvars, int nnonz, int blocksize, int *nvarnonz, int **col, int **row, SCIP_Real **val, SCIP_VAR **vars, int constnnonz, int *constcol, int *constrow, SCIP_Real *constval)
Definition: cons_sdp.c:5150
SCIPconsSdpGetMaxSdpCoef
SCIP_Real SCIPconsSdpGetMaxSdpCoef(SCIP *scip, SCIP_CONS *cons)
Definition: cons_sdp.c:4972
scaleRowsMatrix
static SCIP_RETCODE scaleRowsMatrix(int blocksize, SCIP_Real *matrix, SCIP_Real *scale)
Definition: cons_sdp.c:196
addTwoMinorLinConstraints
static SCIP_RETCODE addTwoMinorLinConstraints(SCIP *scip, SCIP_CONS **conss, int nconss, int *naddconss)
Definition: cons_sdp.c:1117
SCIPconsSdpCompLowerTriangPos
int SCIPconsSdpCompLowerTriangPos(int i, int j)
Definition: cons_sdp.c:4585
SCIPconsSdpGetMaxEVSubmat
SCIP_RETCODE SCIPconsSdpGetMaxEVSubmat(SCIP_CONS *cons, int **maxevsubmat)
Definition: cons_sdp.c:5114
SCIPlapackMatrixVectorMult
SCIP_RETCODE SCIPlapackMatrixVectorMult(int nrows, int ncols, SCIP_Real *matrix, SCIP_Real *vector, SCIP_Real *result)
Definition: lapack_interface.c:354
SCIPlapackComputeIthEigenvalue
SCIP_EXPORT SCIP_RETCODE SCIPlapackComputeIthEigenvalue(BMS_BUFMEM *bufmem, SCIP_Bool geteigenvectors, int n, SCIP_Real *A, int i, SCIP_Real *eigenvalue, SCIP_Real *eigenvector)
SCIPconsSdpGetData
SCIP_RETCODE SCIPconsSdpGetData(SCIP *scip, SCIP_CONS *cons, int *nvars, int *nnonz, int *blocksize, int *arraylength, int *nvarnonz, int **col, int **row, SCIP_Real **val, SCIP_VAR **vars, int *constnnonz, int *constcol, int *constrow, SCIP_Real *constval, SCIP_Bool *rankone, int **maxevsubmat, SCIP_Bool *addedquadcons)
Definition: cons_sdp.c:4603
CONSHDLR_CHECKPRIORITY
#define CONSHDLR_CHECKPRIORITY
Definition: cons_sdp.c:86
SCIP_DECL_CONSENFOPS
static SCIP_DECL_CONSENFOPS(consEnfopsSdp)
Definition: cons_sdp.c:3739
SCIPconsSdpShouldBeRankOne
SCIP_Bool SCIPconsSdpShouldBeRankOne(SCIP_CONS *cons)
Definition: cons_sdp.c:5099
DEFAULT_DIAGZEROIMPLCUTS
#define DEFAULT_DIAGZEROIMPLCUTS
Definition: cons_sdp.c:99
SCIPconsSdpCheckSdpCons
SCIP_RETCODE SCIPconsSdpCheckSdpCons(SCIP *scip, SCIP_CONS *cons, SCIP_SOL *sol, SCIP_Bool printreason, SCIP_RESULT *result)
Definition: cons_sdp.c:540
SdpVarfixer.h
adds the main functionality to fix/unfix/(multi-)aggregate variables by merging two three-tuple-array...
isMatrixRankOne
static SCIP_RETCODE isMatrixRankOne(SCIP *scip, SCIP_CONS *cons, SCIP_SOL *sol, SCIP_Bool *result)
Definition: cons_sdp.c:298
CONSHDLR_ENFOPRIORITY
#define CONSHDLR_ENFOPRIORITY
Definition: cons_sdp.c:85
SCIP_DECL_CONSFREE
static SCIP_DECL_CONSFREE(consFreeSdp)
Definition: cons_sdp.c:3878
SCIP_DECL_CONSPRESOL
static SCIP_DECL_CONSPRESOL(consPresolSdp)
Definition: cons_sdp.c:3124
SCIPincludeConshdlrSdpRank1
SCIP_RETCODE SCIPincludeConshdlrSdpRank1(SCIP *scip)
Definition: cons_sdp.c:4513
SCIP_DECL_CONSEXITPRE
static SCIP_DECL_CONSEXITPRE(consExitpreSdp)
Definition: cons_sdp.c:2867
SCIP_DECL_CONSPRINT
static SCIP_DECL_CONSPRINT(consPrintSdp)
Definition: cons_sdp.c:3999
SCIP_DECL_CONSHDLRCOPY
static SCIP_DECL_CONSHDLRCOPY(conshdlrCopySdp)
Definition: cons_sdp.c:3895
updateVarLocks
static SCIP_RETCODE updateVarLocks(SCIP *scip, SCIP_CONS *cons, int v)
Definition: cons_sdp.c:1599
CONSHDLR_MAXPREROUNDS
#define CONSHDLR_MAXPREROUNDS
Definition: cons_sdp.c:91
CONSHDLRRANK1_NAME
#define CONSHDLRRANK1_NAME
Definition: cons_sdp.c:81
DEFAULT_UPGRADQUADCONSS
#define DEFAULT_UPGRADQUADCONSS
Definition: cons_sdp.c:103
SCIP_DECL_CONSINITPRE
static SCIP_DECL_CONSINITPRE(consInitpreSdp)
Definition: cons_sdp.c:2717
SCIP_DECL_CONSINITSOL
static SCIP_DECL_CONSINITSOL(consInitsolSdp)
Definition: cons_sdp.c:2917
SCIP_DECL_CONSPARSE
static SCIP_DECL_CONSPARSE(consParseSdp)
Definition: cons_sdp.c:4140
CONSHDLRRANK1_DESC
#define CONSHDLRRANK1_DESC
Definition: cons_sdp.c:82
diagZeroImpl
static SCIP_RETCODE diagZeroImpl(SCIP *scip, SCIP_CONS **conss, int nconss, int *naddconss)
Definition: cons_sdp.c:869
PARSE_STARTSIZE
#define PARSE_STARTSIZE
Definition: cons_sdp.c:96
SCIPconsSdpGetBlocksize
int SCIPconsSdpGetBlocksize(SCIP *scip, SCIP_CONS *cons)
Definition: cons_sdp.c:4735
DEFAULT_UPGRADEKEEPQUAD
#define DEFAULT_UPGRADEKEEPQUAD
Definition: cons_sdp.c:104
SCIP_DECL_CONSSEPASOL
static SCIP_DECL_CONSSEPASOL(consSepasolSdp)
Definition: cons_sdp.c:3803
SCIPsdpVarfixerMergeArrays
SCIP_RETCODE SCIPsdpVarfixerMergeArrays(BMS_BLKMEM *blkmem, SCIP_Real epsilon, int *originrow, int *origincol, SCIP_Real *originval, int originlength, SCIP_Bool originsorted, SCIP_Real scalar, int *targetrow, int *targetcol, SCIP_Real *targetval, int *targetlength, int targetmemory)
Definition: SdpVarfixer.c:87
diagGEzero
static SCIP_RETCODE diagGEzero(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_CONS **conss, int nconss, int *naddconss, int *nchgbds, SCIP_RESULT *result)
Definition: cons_sdp.c:698
SCIP_DECL_CONSGETVARS
static SCIP_DECL_CONSGETVARS(consGetVarsSdp)
Definition: cons_sdp.c:4367
CONSHDLR_SEPAPRIORITY
#define CONSHDLR_SEPAPRIORITY
Definition: cons_sdp.c:84
SCIPlapackLinearSolve
SCIP_RETCODE SCIPlapackLinearSolve(BMS_BUFMEM *bufmem, int m, int n, SCIP_Real *A, SCIP_Real *b, SCIP_Real *x)
Definition: lapack_interface.c:438
SCIPconsSdpGetNNonz
SCIP_RETCODE SCIPconsSdpGetNNonz(SCIP *scip, SCIP_CONS *cons, int *nnonz, int *constnnonz)
Definition: cons_sdp.c:4710
lapack_interface.h
interface methods for eigenvector computation and matrix multiplication using openblas ...
fixAndAggrVars
static SCIP_RETCODE fixAndAggrVars(SCIP *scip, SCIP_CONS **conss, int nconss, SCIP_Bool aggregate)
Definition: cons_sdp.c:1966
SCIP_DECL_CONSLOCK
static SCIP_DECL_CONSLOCK(consLockSdp)
Definition: cons_sdp.c:2735
SCIPlapackComputeEigenvectorDecomposition
SCIP_RETCODE SCIPlapackComputeEigenvectorDecomposition(BMS_BUFMEM *bufmem, int n, SCIP_Real *A, SCIP_Real *eigenvalues, SCIP_Real *eigenvectors)
Definition: lapack_interface.c:250
name
char name[SCIP_MAXSTRLEN]
Definition: sdpisolver_mosek.c:267
SCIP_DECL_CONSGETNVARS
static SCIP_DECL_CONSGETNVARS(consGetNVarsSdp)
Definition: cons_sdp.c:4401
SCIP_DECL_CONSCOPY
static SCIP_DECL_CONSCOPY(consCopySdp)
Definition: cons_sdp.c:3925
DEFAULT_TWOMINORLINCONSS
#define DEFAULT_TWOMINORLINCONSS
Definition: cons_sdp.c:100
SCIPconsSdpGuessInitialPoint
SCIP_RETCODE SCIPconsSdpGuessInitialPoint(SCIP *scip, SCIP_CONS *cons, SCIP_Real *lambdastar)
Definition: cons_sdp.c:4860
SCIP_DECL_QUADCONSUPGD
static SCIP_DECL_QUADCONSUPGD(consQuadConsUpgdSdp)
Definition: cons_sdp.c:2304
multiaggrVar
static SCIP_RETCODE multiaggrVar(SCIP *scip, SCIP_CONS *cons, int v, SCIP_VAR **aggrvars, SCIP_Real *scalars, int naggrvars, SCIP_Real constant, int *savedcol, int *savedrow, SCIP_Real *savedval, int *nfixednonz, int *vararraylength)
Definition: cons_sdp.c:1747
unlockVar
static SCIP_RETCODE unlockVar(SCIP *scip, SCIP_CONSDATA *consdata, int v)
Definition: cons_sdp.c:1566
cutUsingEigenvector
static SCIP_RETCODE cutUsingEigenvector(SCIP *scip, SCIP_CONS *cons, SCIP_SOL *sol, SCIP_Real *coeff, SCIP_Real *lhs)
Definition: cons_sdp.c:471
SCIP_DECL_CONSENFOLP
static SCIP_DECL_CONSENFOLP(consEnfolpSdp)
Definition: cons_sdp.c:3778
Generated by   doxygen 1.8.11