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-2016 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 (version number) */
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  assert(N > 0);
229  assert(A != NULL);
230  assert(b != NULL);
231  assert(x != NULL);
232  assert(success != NULL);
233 
234  /* call SCIPsolveLinearProb3() for performance reasons */
235  if( N == 3 )
236  {
237  SCIP_CALL( SCIPsolveLinearProb3(A, b, x, success) );
238  return SCIP_OKAY;
239  }
240 
241  *success = TRUE;
242 
243  LU = NULL;
244  y = NULL;
245  pivot = NULL;
246 
247  /* copy arrays */
248  SCIP_ALLOC( BMSduplicateMemoryArray(&LU, A, N*N) ); /*lint !e647*/
249  SCIP_ALLOC( BMSduplicateMemoryArray(&y, b, N) );
250  SCIP_ALLOC( BMSallocMemoryArray(&pivot, N) );
251 
252  /* initialize values */
253  for( k = 0; k < N; ++k )
254  pivot[k] = k;
255 
256  /* first step: compute LU factorization */
257  for( k = 0; k < N; ++k )
258  {
259  int p;
260  int i;
261 
262  p = k;
263  for( i = k+1; i < N; ++i )
264  {
265  if( ABS(LU[ ENTRY(pivot[i],k) ]) > ABS( LU[ ENTRY(pivot[p],k) ]) )
266  p = i;
267  }
268 
269  if( ABS(LU[ ENTRY(pivot[p],k) ]) < 1e-08 )
270  {
271  SCIPerrorMessage("Error in nlpi_ipopt_dummy - matrix is singular!\n");
272  *success = FALSE;
273  goto TERMINATE;
274  }
275 
276  if( p != k )
277  {
278  int tmp;
279 
280  tmp = pivot[k];
281  pivot[k] = pivot[p];
282  pivot[p] = tmp;
283  }
284 
285  for( i = k+1; i < N; ++i )
286  {
287  SCIP_Real m;
288  int j;
289 
290  m = LU[ ENTRY(pivot[i],k) ] / LU[ ENTRY(pivot[k],k) ];
291 
292  for( j = k+1; j < N; ++j )
293  LU[ ENTRY(pivot[i],j) ] -= m * LU[ ENTRY(pivot[k],j) ];
294 
295  LU[ ENTRY(pivot[i],k) ] = m;
296  }
297  }
298 
299  /* second step: forward substitution */
300  y[0] = b[pivot[0]];
301 
302  for( k = 1; k < N; ++k )
303  {
304  SCIP_Real s;
305  int j;
306 
307  s = b[pivot[k]];
308  for( j = 0; j < k; ++j )
309  {
310  s -= LU[ ENTRY(pivot[k],j) ] * y[j];
311  }
312  y[k] = s;
313  }
314 
315  /* third step: backward substitution */
316  x[N-1] = y[N-1] / LU[ ENTRY(pivot[N-1],N-1) ];
317  for( k = N-2; k >= 0; --k )
318  {
319  SCIP_Real s;
320  int j;
321 
322  s = y[k];
323  for( j = k+1; j < N; ++j )
324  {
325  s -= LU[ ENTRY(pivot[k],j) ] * x[j];
326  }
327  x[k] = s / LU[ ENTRY(pivot[k],k) ];
328  }
329 
330  TERMINATE:
331  /* free arrays */
332  BMSfreeMemoryArray(&pivot);
333  BMSfreeMemoryArray(&y);
334  BMSfreeMemoryArray(&LU);
335 
336  return SCIP_OKAY;
337 }
#define BMSfreeMemoryArray(ptr)
Definition: memory.h:102
SCIP_RETCODE SCIPcreateNlpSolverIpopt(BMS_BLKMEM *blkmem, SCIP_NLPI **nlpi)
void * SCIPgetIpoptApplicationPointerIpopt(SCIP_NLPIPROBLEM *nlpiproblem)
SCIP_RETCODE LapackDsyev(SCIP_Bool computeeigenvectors, int N, SCIP_Real *a, SCIP_Real *w)
#define ENTRY(i, j)
#define NULL
Definition: lpi_spx.cpp:130
void SCIPsetModifiedDefaultSettingsIpopt(SCIP_NLPI *nlpi, const char *optionsstring)
#define FALSE
Definition: def.h:56
#define TRUE
Definition: def.h:55
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:53
#define SCIP_CALL(x)
Definition: def.h:266
SCIP_RETCODE SCIPsolveLinearProb(int N, SCIP_Real *A, SCIP_Real *b, SCIP_Real *x, SCIP_Bool *success)
#define SCIPerrorMessage
Definition: pub_message.h:45
void * SCIPgetNlpiOracleIpopt(SCIP_NLPIPROBLEM *nlpiproblem)
#define BMSduplicateMemoryArray(ptr, source, num)
Definition: memory.h:98
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:53
#define BMScopyMemoryArray(ptr, source, num)
Definition: memory.h:89
public methods for message output
SCIP_Bool SCIPisIpoptAvailableIpopt(void)
#define SCIP_Real
Definition: def.h:127
const char * SCIPgetSolverDescIpopt(void)
#define BMSallocMemoryArray(ptr, num)
Definition: memory.h:78
struct BMS_BlkMem BMS_BLKMEM
Definition: memory.h:392
const char * SCIPgetSolverNameIpopt(void)
#define SCIP_ALLOC(x)
Definition: def.h:277
#define SCIPABORT()
Definition: def.h:238