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-2018 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not email to scip@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file nlpi_ipopt_dummy.c
17  * @brief dummy Ipopt NLP interface for the case that Ipopt is not available
18  * @author Stefan Vigerske
19  * @author Benjamin Müller
20  *
21  * This code has been separate from nlpi_ipopt.cpp, so the SCIP build system recognizes it as pure C code,
22  * thus the linker does not need to be changed to C++.
23  */
24 
25 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
26 
27 #include "scip/pub_message.h"
28 #include "nlpi/nlpi_ipopt.h"
29 
30 /** create solver interface for Ipopt solver */
32  BMS_BLKMEM* blkmem, /**< block memory data structure */
33  SCIP_NLPI** nlpi /**< pointer to buffer for nlpi address */
34  )
35 {
36  assert(nlpi != NULL);
37 
38  *nlpi = NULL;
39 
40  return SCIP_OKAY;
41 } /*lint !e715*/
42 
43 /** gets string that identifies Ipopt (version number) */
44 const char* SCIPgetSolverNameIpopt(void)
45 {
46  return "";
47 }
48 
49 /** gets string that describes Ipopt */
50 const char* SCIPgetSolverDescIpopt(void)
51 {
52  return "";
53 }
54 
55 /** returns whether Ipopt is available, i.e., whether it has been linked in */
57 {
58  return FALSE;
59 }
60 
61 /** gives a pointer to the IpoptApplication object stored in Ipopt-NLPI's NLPI problem data structure */
63  SCIP_NLPIPROBLEM* nlpiproblem /**< NLP problem of Ipopt-NLPI */
64  )
65 {
66  SCIPerrorMessage("Ipopt not available!\n");
67  SCIPABORT();
68  return NULL; /*lint !e527*/
69 } /*lint !e715*/
70 
71 /** gives a pointer to the NLPIORACLE object stored in Ipopt-NLPI's NLPI problem data structure */
73  SCIP_NLPIPROBLEM* nlpiproblem /**< NLP problem of Ipopt-NLPI */
74  )
75 {
76  SCIPerrorMessage("Ipopt not available!\n");
77  SCIPABORT();
78  return NULL; /*lint !e527*/
79 } /*lint !e715*/
80 
81 /** sets modified default settings that are used when setting up an Ipopt problem
82  *
83  * Do not forget to add a newline after the last option in optionsstring.
84  */
86  SCIP_NLPI* nlpi, /**< Ipopt NLP interface */
87  const char* optionsstring /**< string with options as in Ipopt options file */
88  )
89 {
90  SCIPerrorMessage("Ipopt not available!\n");
91  SCIPABORT();
92 } /*lint !e715*/
93 
94 /** Calls Lapacks Dsyev routine to compute eigenvalues and eigenvectors of a dense matrix.
95  * It's here, because Ipopt is linked against Lapack.
96  */
98  SCIP_Bool computeeigenvectors,/**< should also eigenvectors should be computed ? */
99  int N, /**< dimension */
100  SCIP_Real* a, /**< matrix data on input (size N*N); eigenvectors on output if computeeigenvectors == TRUE */
101  SCIP_Real* w /**< buffer to store eigenvalues (size N) */
102  )
103 {
104  SCIPerrorMessage("Ipopt not available, cannot use it's Lapack link!\n");
105  return SCIP_ERROR;
106 } /*lint !e715*/
107 
108 /* easier access to the entries of A */
109 #define ENTRY(i,j) (N * (j) + (i))
110 
111 /* solves a linear problem of the form Ax = b for a regular 3*3 matrix A */
112 static
114  SCIP_Real* A, /**< matrix data on input (size 3*3); filled column-wise */
115  SCIP_Real* b, /**< right hand side vector (size 3) */
116  SCIP_Real* x, /**< buffer to store solution (size 3) */
117  SCIP_Bool* success /**< pointer to store if the solving routine was successful */
118  )
119 {
120  SCIP_Real LU[9];
121  SCIP_Real y[3];
122  int pivot[3] = {0, 1, 2};
123  const int N = 3;
124  int k;
125 
126  assert(A != NULL);
127  assert(b != NULL);
128  assert(x != NULL);
129  assert(success != NULL);
130 
131  *success = TRUE;
132 
133  /* copy arrays */
134  BMScopyMemoryArray(LU, A, N*N);
135  BMScopyMemoryArray(y, b, N);
136 
137  /* first step: compute LU factorization */
138  for( k = 0; k < N; ++k )
139  {
140  int p;
141  int i;
142 
143  p = k;
144  for( i = k+1; i < N; ++i )
145  {
146  if( ABS(LU[ ENTRY(pivot[i],k) ]) > ABS( LU[ ENTRY(pivot[p],k) ]) )
147  p = i;
148  }
149 
150  if( ABS(LU[ ENTRY(pivot[p],k) ]) < 1e-08 )
151  {
152  SCIPerrorMessage("Error in nlpi_ipopt_dummy - matrix is singular!\n");
153  *success = FALSE;
154  return SCIP_OKAY;
155  }
156 
157  if( p != k )
158  {
159  int tmp;
160 
161  tmp = pivot[k];
162  pivot[k] = pivot[p];
163  pivot[p] = tmp;
164  }
165 
166  for( i = k+1; i < N; ++i )
167  {
168  SCIP_Real m;
169  int j;
170 
171  m = LU[ ENTRY(pivot[i],k) ] / LU[ ENTRY(pivot[k],k) ];
172 
173  for( j = k+1; j < N; ++j )
174  LU[ ENTRY(pivot[i],j) ] -= m * LU[ ENTRY(pivot[k],j) ];
175 
176  LU[ ENTRY(pivot[i],k) ] = m;
177  }
178  }
179 
180  /* second step: forward substitution */
181  y[0] = b[pivot[0]];
182 
183  for( k = 1; k < N; ++k )
184  {
185  SCIP_Real s;
186  int j;
187 
188  s = b[pivot[k]];
189  for( j = 0; j < k; ++j )
190  {
191  s -= LU[ ENTRY(pivot[k],j) ] * y[j];
192  }
193  y[k] = s;
194  }
195 
196  /* third step: backward substitution */
197  x[N-1] = y[N-1] / LU[ ENTRY(pivot[N-1],N-1) ];
198  for( k = N-2; k >= 0; --k )
199  {
200  SCIP_Real s;
201  int j;
202 
203  s = y[k];
204  for( j = k+1; j < N; ++j )
205  {
206  s -= LU[ ENTRY(pivot[k],j) ] * x[j];
207  }
208  x[k] = s / LU[ ENTRY(pivot[k],k) ];
209  }
210 
211  return SCIP_OKAY;
212 }
213 
214 /* solves a linear problem of the form Ax = b for a regular matrix A */
216  int N, /**< dimension */
217  SCIP_Real* A, /**< matrix data on input (size N*N); filled column-wise */
218  SCIP_Real* b, /**< right hand side vector (size N) */
219  SCIP_Real* x, /**< buffer to store solution (size N) */
220  SCIP_Bool* success /**< pointer to store if the solving routine was successful */
221  )
222 {
223  SCIP_Real* LU;
224  SCIP_Real* y;
225  int* pivot;
226  int k;
227 
228  SCIP_RETCODE retcode = SCIP_OKAY;
229 
230  assert(N > 0);
231  assert(A != NULL);
232  assert(b != NULL);
233  assert(x != NULL);
234  assert(success != NULL);
235 
236  /* call SCIPsolveLinearProb3() for performance reasons */
237  if( N == 3 )
238  {
239  SCIP_CALL( SCIPsolveLinearProb3(A, b, x, success) );
240  return SCIP_OKAY;
241  }
242 
243  *success = TRUE;
244 
245  LU = NULL;
246  y = NULL;
247  pivot = NULL;
248 
249  /* copy arrays */
250  SCIP_ALLOC_TERMINATE( retcode, BMSduplicateMemoryArray(&LU, A, N*N), TERMINATE ); /*lint !e647*/
251  SCIP_ALLOC_TERMINATE( retcode, BMSduplicateMemoryArray(&y, b, N), TERMINATE );
252  SCIP_ALLOC_TERMINATE( retcode, BMSallocMemoryArray(&pivot, N), TERMINATE );
253 
254  /* initialize values */
255  for( k = 0; k < N; ++k )
256  pivot[k] = k;
257 
258  /* first step: compute LU factorization */
259  for( k = 0; k < N; ++k )
260  {
261  int p;
262  int i;
263 
264  p = k;
265  for( i = k+1; i < N; ++i )
266  {
267  if( ABS(LU[ ENTRY(pivot[i],k) ]) > ABS( LU[ ENTRY(pivot[p],k) ]) )
268  p = i;
269  }
270 
271  if( ABS(LU[ ENTRY(pivot[p],k) ]) < 1e-08 )
272  {
273  SCIPerrorMessage("Error in nlpi_ipopt_dummy - matrix is singular!\n");
274  *success = FALSE;
275  goto TERMINATE;
276  }
277 
278  if( p != k )
279  {
280  int tmp;
281 
282  tmp = pivot[k];
283  pivot[k] = pivot[p];
284  pivot[p] = tmp;
285  }
286 
287  for( i = k+1; i < N; ++i )
288  {
289  SCIP_Real m;
290  int j;
291 
292  m = LU[ ENTRY(pivot[i],k) ] / LU[ ENTRY(pivot[k],k) ];
293 
294  for( j = k+1; j < N; ++j )
295  LU[ ENTRY(pivot[i],j) ] -= m * LU[ ENTRY(pivot[k],j) ];
296 
297  LU[ ENTRY(pivot[i],k) ] = m;
298  }
299  }
300 
301  /* second step: forward substitution */
302  y[0] = b[pivot[0]];
303 
304  for( k = 1; k < N; ++k )
305  {
306  SCIP_Real s;
307  int j;
308 
309  s = b[pivot[k]];
310  for( j = 0; j < k; ++j )
311  {
312  s -= LU[ ENTRY(pivot[k],j) ] * y[j];
313  }
314  y[k] = s;
315  }
316 
317  /* third step: backward substitution */
318  x[N-1] = y[N-1] / LU[ ENTRY(pivot[N-1],N-1) ];
319  for( k = N-2; k >= 0; --k )
320  {
321  SCIP_Real s;
322  int j;
323 
324  s = y[k];
325  for( j = k+1; j < N; ++j )
326  {
327  s -= LU[ ENTRY(pivot[k],j) ] * x[j];
328  }
329  x[k] = s / LU[ ENTRY(pivot[k],k) ];
330  }
331 
332  TERMINATE:
333  /* free arrays */
334  BMSfreeMemoryArrayNull(&pivot);
337 
338  return retcode;
339 }
SCIP_Bool SCIPisIpoptAvailableIpopt(void)
void * SCIPgetNlpiOracleIpopt(SCIP_NLPIPROBLEM *nlpiproblem)
#define ENTRY(i, j)
#define BMSfreeMemoryArrayNull(ptr)
Definition: memory.h:130
void SCIPsetModifiedDefaultSettingsIpopt(SCIP_NLPI *nlpi, const char *optionsstring)
#define FALSE
Definition: def.h:64
#define TRUE
Definition: def.h:63
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:53
#define BMSallocMemoryArray(ptr, num)
Definition: memory.h:105
const char * SCIPgetSolverNameIpopt(void)
#define SCIPerrorMessage
Definition: pub_message.h:45
#define SCIP_CALL(x)
Definition: def.h:350
SCIP_RETCODE LapackDsyev(SCIP_Bool computeeigenvectors, int N, SCIP_Real *a, SCIP_Real *w)
#define BMSduplicateMemoryArray(ptr, source, num)
Definition: memory.h:125
Ipopt NLP interface.
static SCIP_RETCODE SCIPsolveLinearProb3(SCIP_Real *A, SCIP_Real *b, SCIP_Real *x, SCIP_Bool *success)
#define SCIP_Bool
Definition: def.h:61
#define BMScopyMemoryArray(ptr, source, num)
Definition: memory.h:116
SCIP_RETCODE SCIPcreateNlpSolverIpopt(BMS_BLKMEM *blkmem, SCIP_NLPI **nlpi)
SCIP_RETCODE SCIPsolveLinearProb(int N, SCIP_Real *A, SCIP_Real *b, SCIP_Real *x, SCIP_Bool *success)
void * SCIPgetIpoptApplicationPointerIpopt(SCIP_NLPIPROBLEM *nlpiproblem)
public methods for message output
const char * SCIPgetSolverDescIpopt(void)
#define SCIP_Real
Definition: def.h:149
struct BMS_BlkMem BMS_BLKMEM
Definition: memory.h:419
#define SCIP_ALLOC_TERMINATE(retcode, x, TERM)
Definition: def.h:381
#define SCIPABORT()
Definition: def.h:322