I have a quadratic optimization problem with linear constraints that I want to solve using SCIP. The optimization matrix that I want to be minimized is positive semi-definite (it is the variance of certain variables, to be precise). I have the problem in a file in CPLEX LP format and when I optimize in SCIP, I get the message
Quadratic constraint handler does not have LAPACK for eigenvalue computation. Will assume
that matrices (with size > 2x2) are indefinite.
So SCIP starts optimization assuming that the matrix is indefinite and takes a large amount of time. I have installed LAPACK and even copied liblapack.a file in the lib folder where the SCIP source and binaries are and reinstalled SCIP. But, I keep getting the above message.
Is there a way to make SCIP use the LAPACK library? I believe the optimization will be really fast, if SCIP can figure out that the matrix is positive semi-definite.
If you feel like patching up SCIP a bit to use your Lapack lib without providing a full Ipopt (though it's relatively easy to build on *nix and could help performance, as mattmilten pointed out), here is a patch that you could try out:
diff --git a/src/scip/cons_quadratic.c b/src/scip/cons_quadratic.c
index 93ba359..795bade 100644
--- a/src/scip/cons_quadratic.c
+++ b/src/scip/cons_quadratic.c
@@ -46,7 +46,7 @@
#include "scip/heur_trysol.h"
#include "scip/debug.h"
#include "nlpi/nlpi.h"
-#include "nlpi/nlpi_ipopt.h"
+/*#include "nlpi/nlpi_ipopt.h" */
/* constraint handler properties */
#define CONSHDLR_NAME "quadratic"
@@ -4257,6 +4257,71 @@ void checkCurvatureEasy(
*determined = FALSE;
}
+#define F77_FUNC(a,A) a##_
+
+ /** LAPACK Fortran subroutine DSYEV */
+ void F77_FUNC(dsyev,DSYEV)(
+ char* jobz, /**< 'N' to compute eigenvalues only, 'V' to compute eigenvalues and eigenvectors */
+ char* uplo, /**< 'U' if upper triangle of A is stored, 'L' if lower triangle of A is stored */
+ int* n, /**< dimension */
+ double* A, /**< matrix A on entry; orthonormal eigenvectors on exit, if jobz == 'V' and info == 0; if jobz == 'N', then the matrix data is destroyed */
+ int* ldA, /**< leading dimension, probably equal to n */
+ double* W, /**< buffer for the eigenvalues in ascending order */
+ double* WORK, /**< workspace array */
+ int* LWORK, /**< length of WORK; if LWORK = -1, then the optimal workspace size is calculated and returned in WORK(1) */
+ int* info /**< == 0: successful exit; < 0: illegal argument at given position; > 0: failed to converge */
+ );
+
+/** Calls Lapacks Dsyev routine to compute eigenvalues and eigenvectors of a dense matrix.
+ */
+static
+SCIP_RETCODE LapackDsyev(
+ SCIP_Bool computeeigenvectors,/**< should also eigenvectors should be computed ? */
+ int N, /**< dimension */
+ SCIP_Real* a, /**< matrix data on input (size N*N); eigenvectors on output if computeeigenvectors == TRUE */
+ SCIP_Real* w /**< buffer to store eigenvalues (size N) */
+ )
+{
+ int INFO;
+ char JOBZ = computeeigenvectors ? 'V' : 'N';
+ char UPLO = 'L';
+ int LDA = N;
+ double* WORK = NULL;
+ int LWORK;
+ double WORK_PROBE;
+ int i;
+
+ /* First we find out how large LWORK should be */
+ LWORK = -1;
+ F77_FUNC(dsyev,DSYEV)(&JOBZ, &UPLO, &N, a, &LDA, w, &WORK_PROBE, &LWORK, &INFO);
+ if( INFO != 0 )
+ {
+ SCIPerrorMessage("There was an error when calling DSYEV. INFO = %d\n", INFO);
+ return SCIP_ERROR;
+ }
+
+ LWORK = (int) WORK_PROBE;
+ assert(LWORK > 0);
+
+ SCIP_ALLOC( BMSallocMemoryArray(&WORK, LWORK) );
+
+ for( i = 0; i < LWORK; ++i )
+ WORK[i] = i;
+
+ F77_FUNC(dsyev,DSYEV)(&JOBZ, &UPLO, &N, a, &LDA, w, WORK, &LWORK, &INFO);
+
+ BMSfreeMemoryArray(&WORK);
+
+ if( INFO != 0 )
+ {
+ SCIPerrorMessage("There was an error when calling DSYEV. INFO = %d\n", INFO);
+ return SCIP_ERROR;
+ }
+
+ return SCIP_OKAY;
+}
+
+
/** checks a quadratic constraint for convexity and/or concavity */
static
SCIP_RETCODE checkCurvature(
@@ -4343,7 +4408,7 @@ SCIP_RETCODE checkCurvature(
return SCIP_OKAY;
}
- if( SCIPisIpoptAvailableIpopt() )
+ if( TRUE )
{
for( i = 0; i < consdata->nbilinterms; ++i )
{
@@ -4479,7 +4544,7 @@ SCIP_RETCODE checkFactorable(
return SCIP_OKAY;
/* need routine to compute eigenvalues/eigenvectors */
- if( !SCIPisIpoptAvailableIpopt() )
+ if( !TRUE )
return SCIP_OKAY;
SCIP_CALL( consdataSortQuadVarTerms(scip, consdata) );
@@ -9395,7 +9460,7 @@ SCIP_DECL_CONSINITSOL(consInitsolQuadratic)
SCIP_CALL( SCIPcatchEvent(scip, SCIP_EVENTTYPE_SOLFOUND, eventhdlr, (SCIP_EVENTDATA*)conshdlr, &conshdlrdata->newsoleventfilterpos) );
}
- if( nconss != 0 && !SCIPisIpoptAvailableIpopt() && !SCIPisInRestart(scip) )
+ if( nconss != 0 && !TRUE && !SCIPisInRestart(scip) )
{
SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, "Quadratic constraint handler does not have LAPACK for eigenvalue computation. Will assume that matrices (with size > 2x2) are indefinite.\n");
}
Use USRLDFLAGS="-llapack -lblas" with make.