Scippy

SCIP

Solving Constraint Integer Programs

nlpi_ipopt_dummy.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB) */
7 /* */
8 /* Licensed under the Apache License, Version 2.0 (the "License"); */
9 /* you may not use this file except in compliance with the License. */
10 /* You may obtain a copy of the License at */
11 /* */
12 /* http://www.apache.org/licenses/LICENSE-2.0 */
13 /* */
14 /* Unless required by applicable law or agreed to in writing, software */
15 /* distributed under the License is distributed on an "AS IS" BASIS, */
16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17 /* See the License for the specific language governing permissions and */
18 /* limitations under the License. */
19 /* */
20 /* You should have received a copy of the Apache-2.0 license */
21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */
22 /* */
23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24 
25 /**@file nlpi_ipopt_dummy.c
26  * @ingroup DEFPLUGINS_NLPI
27  * @brief dummy Ipopt NLP interface for the case that Ipopt is not available
28  * @author Stefan Vigerske
29  * @author Benjamin Müller
30  *
31  * This code has been separated from nlpi_ipopt.cpp, so the SCIP build system recognizes it as pure C code,
32  * thus the linker does not need to be changed to C++.
33  */
34 
35 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
36 
37 #include "scip/pub_message.h"
38 #include "scip/nlpi_ipopt.h"
39 #include "blockmemshell/memory.h"
40 
41 /** create solver interface for Ipopt solver and includes it into SCIP, if Ipopt is available */ /*lint -e{715}*/
43  SCIP* scip /**< SCIP data structure */
44  )
45 {
46  return SCIP_OKAY;
47 } /*lint !e715*/
48 
49 /** gets string that identifies Ipopt (version number) */
50 const char* SCIPgetSolverNameIpopt(void)
51 {
52  return "";
53 }
54 
55 /** gets string that describes Ipopt */
56 const char* SCIPgetSolverDescIpopt(void)
57 {
58  return "";
59 }
60 
61 /** returns whether Ipopt is available, i.e., whether it has been linked in */
63 {
64  return FALSE;
65 }
66 
67 /** gives a pointer to the NLPIORACLE object stored in Ipopt-NLPI's NLPI problem data structure */ /*lint -e715*/
69  SCIP_NLPIPROBLEM* nlpiproblem /**< NLP problem of Ipopt-NLPI */
70  )
71 {
72  SCIPerrorMessage("Ipopt not available!\n");
73  SCIPABORT();
74  return NULL; /*lint !e527*/
75 } /*lint !e715*/
76 
77 /** Calls Lapacks Dsyev routine to compute eigenvalues and eigenvectors of a dense matrix.
78  * It's here, because Ipopt is linked against Lapack.
79  */ /*lint -e715*/
81  SCIP_Bool computeeigenvectors,/**< should also eigenvectors should be computed ? */
82  int N, /**< dimension */
83  SCIP_Real* a, /**< matrix data on input (size N*N); eigenvectors on output if computeeigenvectors == TRUE */
84  SCIP_Real* w /**< buffer to store eigenvalues (size N) */
85  )
86 {
87  SCIPerrorMessage("Ipopt not available, cannot use it's Lapack link!\n");
88  return SCIP_ERROR;
89 } /*lint !e715*/
90 
91 /* easier access to the entries of A */
92 #define ENTRY(i,j) (N * (j) + (i))
93 
94 /* solves a linear problem of the form Ax = b for a regular 3*3 matrix A */
95 static
97  SCIP_Real* A, /**< matrix data on input (size 3*3); filled column-wise */
98  SCIP_Real* b, /**< right hand side vector (size 3) */
99  SCIP_Real* x, /**< buffer to store solution (size 3) */
100  SCIP_Bool* success /**< pointer to store if the solving routine was successful */
101  )
102 {
103  SCIP_Real LU[9];
104  SCIP_Real y[3];
105  int pivot[3] = {0, 1, 2};
106  const int N = 3;
107  int k;
108 
109  assert(A != NULL);
110  assert(b != NULL);
111  assert(x != NULL);
112  assert(success != NULL);
113 
114  *success = TRUE;
115 
116  /* copy arrays */
117  BMScopyMemoryArray(LU, A, N*N);
118  BMScopyMemoryArray(y, b, N);
119 
120  /* first step: compute LU factorization */
121  for( k = 0; k < N; ++k )
122  {
123  int p;
124  int i;
125 
126  p = k;
127  for( i = k+1; i < N; ++i )
128  {
129  if( ABS(LU[ ENTRY(pivot[i],k) ]) > ABS( LU[ ENTRY(pivot[p],k) ]) )
130  p = i;
131  }
132 
133  if( ABS(LU[ ENTRY(pivot[p],k) ]) < 1e-08 )
134  {
135  /* SCIPerrorMessage("Error in nlpi_ipopt_dummy - matrix is singular!\n"); */
136  *success = FALSE;
137  return SCIP_OKAY;
138  }
139 
140  if( p != k )
141  {
142  int tmp;
143 
144  tmp = pivot[k];
145  pivot[k] = pivot[p];
146  pivot[p] = tmp;
147  }
148 
149  for( i = k+1; i < N; ++i )
150  {
151  SCIP_Real m;
152  int j;
153 
154  m = LU[ ENTRY(pivot[i],k) ] / LU[ ENTRY(pivot[k],k) ];
155 
156  for( j = k+1; j < N; ++j )
157  LU[ ENTRY(pivot[i],j) ] -= m * LU[ ENTRY(pivot[k],j) ];
158 
159  LU[ ENTRY(pivot[i],k) ] = m;
160  }
161  }
162 
163  /* second step: forward substitution */
164  y[0] = b[pivot[0]];
165 
166  for( k = 1; k < N; ++k )
167  {
168  SCIP_Real s;
169  int j;
170 
171  s = b[pivot[k]];
172  for( j = 0; j < k; ++j )
173  {
174  s -= LU[ ENTRY(pivot[k],j) ] * y[j];
175  }
176  y[k] = s;
177  }
178 
179  /* third step: backward substitution */
180  x[N-1] = y[N-1] / LU[ ENTRY(pivot[N-1],N-1) ];
181  for( k = N-2; k >= 0; --k )
182  {
183  SCIP_Real s;
184  int j;
185 
186  s = y[k];
187  for( j = k+1; j < N; ++j )
188  {
189  s -= LU[ ENTRY(pivot[k],j) ] * x[j];
190  }
191  x[k] = s / LU[ ENTRY(pivot[k],k) ];
192  }
193 
194  return SCIP_OKAY;
195 }
196 
197 /* solves a linear problem of the form Ax = b for a regular matrix A */
199  int N, /**< dimension */
200  SCIP_Real* A, /**< matrix data on input (size N*N); filled column-wise */
201  SCIP_Real* b, /**< right hand side vector (size N) */
202  SCIP_Real* x, /**< buffer to store solution (size N) */
203  SCIP_Bool* success /**< pointer to store if the solving routine was successful */
204  )
205 {
206  SCIP_Real* LU;
207  SCIP_Real* y;
208  int* pivot;
209  int k;
210 
211  SCIP_RETCODE retcode = SCIP_OKAY;
212 
213  assert(N > 0);
214  assert(A != NULL);
215  assert(b != NULL);
216  assert(x != NULL);
217  assert(success != NULL);
218 
219  /* call solveLinearProb3() for performance reasons */
220  if( N == 3 )
221  {
222  SCIP_CALL( solveLinearProb3(A, b, x, success) );
223  return SCIP_OKAY;
224  }
225 
226  *success = TRUE;
227 
228  LU = NULL;
229  y = NULL;
230  pivot = NULL;
231 
232  /* copy arrays */
233  SCIP_ALLOC_TERMINATE( retcode, BMSduplicateMemoryArray(&LU, A, N*N), TERMINATE ); /*lint !e647*/
234  SCIP_ALLOC_TERMINATE( retcode, BMSduplicateMemoryArray(&y, b, N), TERMINATE );
235  SCIP_ALLOC_TERMINATE( retcode, BMSallocMemoryArray(&pivot, N), TERMINATE );
236 
237  /* initialize values */
238  for( k = 0; k < N; ++k )
239  pivot[k] = k;
240 
241  /* first step: compute LU factorization */
242  for( k = 0; k < N; ++k )
243  {
244  int p;
245  int i;
246 
247  p = k;
248  for( i = k+1; i < N; ++i )
249  {
250  if( ABS(LU[ ENTRY(pivot[i],k) ]) > ABS( LU[ ENTRY(pivot[p],k) ]) )
251  p = i;
252  }
253 
254  if( ABS(LU[ ENTRY(pivot[p],k) ]) < 1e-08 )
255  {
256  /* SCIPerrorMessage("Error in nlpi_ipopt_dummy - matrix is singular!\n"); */
257  *success = FALSE;
258  goto TERMINATE;
259  }
260 
261  if( p != k )
262  {
263  int tmp;
264 
265  tmp = pivot[k];
266  pivot[k] = pivot[p];
267  pivot[p] = tmp;
268  }
269 
270  for( i = k+1; i < N; ++i )
271  {
272  SCIP_Real m;
273  int j;
274 
275  m = LU[ ENTRY(pivot[i],k) ] / LU[ ENTRY(pivot[k],k) ];
276 
277  for( j = k+1; j < N; ++j )
278  LU[ ENTRY(pivot[i],j) ] -= m * LU[ ENTRY(pivot[k],j) ];
279 
280  LU[ ENTRY(pivot[i],k) ] = m;
281  }
282  }
283 
284  /* second step: forward substitution */
285  y[0] = b[pivot[0]];
286 
287  for( k = 1; k < N; ++k )
288  {
289  SCIP_Real s;
290  int j;
291 
292  s = b[pivot[k]];
293  for( j = 0; j < k; ++j )
294  {
295  s -= LU[ ENTRY(pivot[k],j) ] * y[j];
296  }
297  y[k] = s;
298  }
299 
300  /* third step: backward substitution */
301  x[N-1] = y[N-1] / LU[ ENTRY(pivot[N-1],N-1) ];
302  for( k = N-2; k >= 0; --k )
303  {
304  SCIP_Real s;
305  int j;
306 
307  s = y[k];
308  for( j = k+1; j < N; ++j )
309  {
310  s -= LU[ ENTRY(pivot[k],j) ] * x[j];
311  }
312  x[k] = s / LU[ ENTRY(pivot[k],k) ];
313  }
314 
315  TERMINATE:
316  /* free arrays */
317  BMSfreeMemoryArrayNull(&pivot);
320 
321  return retcode;
322 }
SCIP_Bool SCIPisIpoptAvailableIpopt(void)
void * SCIPgetNlpiOracleIpopt(SCIP_NLPIPROBLEM *nlpiproblem)
#define ENTRY(i, j)
#define BMSfreeMemoryArrayNull(ptr)
Definition: memory.h:150
SCIP_RETCODE SCIPincludeNlpSolverIpopt(SCIP *scip)
SCIP_RETCODE SCIPcallLapackDsyevIpopt(SCIP_Bool computeeigenvectors, int N, SCIP_Real *a, SCIP_Real *w)
#define FALSE
Definition: def.h:96
#define TRUE
Definition: def.h:95
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:63
#define BMSallocMemoryArray(ptr, num)
Definition: memory.h:125
SCIP_VAR ** x
Definition: circlepacking.c:63
SCIP_VAR * w
Definition: circlepacking.c:67
const char * SCIPgetSolverNameIpopt(void)
#define SCIPerrorMessage
Definition: pub_message.h:64
#define NULL
Definition: lpi_spx1.cpp:164
SCIP_RETCODE SCIPsolveLinearEquationsIpopt(int N, SCIP_Real *A, SCIP_Real *b, SCIP_Real *x, SCIP_Bool *success)
#define SCIP_CALL(x)
Definition: def.h:394
#define BMSduplicateMemoryArray(ptr, source, num)
Definition: memory.h:145
Ipopt NLP interface.
#define SCIP_Bool
Definition: def.h:93
static SCIP_RETCODE solveLinearProb3(SCIP_Real *A, SCIP_Real *b, SCIP_Real *x, SCIP_Bool *success)
#define BMScopyMemoryArray(ptr, source, num)
Definition: memory.h:136
SCIP_VAR ** b
Definition: circlepacking.c:65
public methods for message output
SCIP_VAR * a
Definition: circlepacking.c:66
const char * SCIPgetSolverDescIpopt(void)
#define SCIP_Real
Definition: def.h:186
SCIP_VAR ** y
Definition: circlepacking.c:64
#define SCIP_ALLOC_TERMINATE(retcode, x, TERM)
Definition: def.h:425
#define SCIPABORT()
Definition: def.h:366
memory allocation routines