Scippy

SCIP

Solving Constraint Integer Programs

lapack_calls.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-2024 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 lapack_calls.h
26 * @brief interface methods for lapack functions
27 * @author Marc Pfetsch
28 *
29 * This file is used to call the LAPACK routine DSYEVR and DGETRF.
30 *
31 * LAPACK can be built with 32- or 64-bit integers, which is not visible to the outside. This interface tries to work
32 * around this issue. Since the Fortran routines are called by reference, they only get a pointer. We always use 64-bit
33 * integers on input, but reduce the output to 32-bit integers. We assume that all sizes can be represented in 32-bit
34 * integers.
35 */
36
37/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
38
39#include <assert.h>
40
41#include "scip/lapack_calls.h"
42
43#include "scip/def.h"
44#include "scip/pub_message.h" /* for debug and error message */
46#include "scip/type_retcode.h"
47#include "scip/nlpi_ipopt.h"
48
49/* turn off lint warnings for whole file: */
50/*lint --e{788,818}*/
51
52
53#ifdef SCIP_WITH_LAPACK
54/* we use 64 bit integers as the base type */
55typedef int64_t LAPACKINTTYPE;
56
57/** transforms a SCIP_Real (that should be integer, but might be off by some numerical error) to an integer by adding 0.5 and rounding down */
58#define SCIP_RealTOINT(x) ((LAPACKINTTYPE) (x + 0.5))
59
60/*
61 * BLAS/LAPACK Calls
62 */
63
64/**@name BLAS/LAPACK Calls */
65/**@{ */
66
67/** Define to a macro mangling the given C identifier (in lower and upper
68 * case), which must not contain underscores, for linking with Fortran. */
69#ifdef FNAME_LCASE_DECOR
70#define F77_FUNC(name,NAME) name ## _
71#endif
72#ifdef FNAME_UCASE_DECOR
73#define F77_FUNC(name,NAME) NAME ## _
74#endif
75#ifdef FNAME_LCASE_NODECOR
76#define F77_FUNC(name,NAME) name
77#endif
78#ifdef FNAME_UCASE_NODECOR
79#define F77_FUNC(name,NAME) NAME
80#endif
81
82/* use backup ... */
83#ifndef F77_FUNC
84#define F77_FUNC(name,NAME) name ## _
85#endif
86
87
88/** LAPACK Fortran subroutine DSYEVR */
89void F77_FUNC(dsyevr, DSYEVR)(char* JOBZ, char* RANGE, char* UPLO,
90 LAPACKINTTYPE* N, SCIP_Real* A, LAPACKINTTYPE* LDA,
91 SCIP_Real* VL, SCIP_Real* VU,
92 LAPACKINTTYPE* IL, LAPACKINTTYPE* IU,
93 SCIP_Real* ABSTOL, LAPACKINTTYPE* M, SCIP_Real* W, SCIP_Real* Z,
94 LAPACKINTTYPE* LDZ, LAPACKINTTYPE* ISUPPZ, SCIP_Real* WORK,
95 LAPACKINTTYPE* LWORK, LAPACKINTTYPE* IWORK, LAPACKINTTYPE* LIWORK,
96 LAPACKINTTYPE* INFO);
97
98/** LAPACK Fortran subroutine DGETRF */
99void F77_FUNC(dgetrf, DGETRF)(LAPACKINTTYPE* M, LAPACKINTTYPE* N, SCIP_Real* A,
100 LAPACKINTTYPE* LDA, LAPACKINTTYPE* IPIV, LAPACKINTTYPE* INFO);
101
102/** LAPACK Fortran subroutine DGETRS */
103void F77_FUNC(dgetrs, DGETRS)(char* TRANS, LAPACKINTTYPE* N, LAPACKINTTYPE* NRHS,
104 SCIP_Real* A, LAPACKINTTYPE* LDA, LAPACKINTTYPE* IPIV, SCIP_Real* B, LAPACKINTTYPE* LDB,
105 LAPACKINTTYPE* INFO);
106
107/** LAPACK Fortran subroutine ivlayer */
108void F77_FUNC(ilaver, ILAVER)(LAPACKINTTYPE* MAJOR, LAPACKINTTYPE* MINOR, LAPACKINTTYPE* PATCH);
109
110/**@} */
111#endif
112
113/*
114 * Functions
115 */
116
117/**@name Functions */
118/**@{ */
119
120/** returns whether Lapack is available, i.e., whether it has been linked in */
122{
124 return TRUE;
125
126#ifdef SCIP_WITH_LAPACK
127 return TRUE;
128#else
129 return FALSE;
130#endif
131}
132
133#ifdef SCIP_WITH_LAPACK
134/** converts a number stored in a int64_t to an int, depending on big- or little endian machines
135 *
136 * We assume that the number actually fits into an int. Thus, if more bits are used, we assume that the number is
137 * negative.
138 */
139static
140int convertToInt(
141 int64_t num /**< number to be converted */
142 )
143{
144 union
145 {
146 int64_t big;
147 int small[2];
148 } work;
149 int checkval = 1;
150
151 assert(sizeof(work) > sizeof(checkval)); /*lint !e506*/
152
153 work.big = num;
154
155 /* if we have a little-endian machine (e.g, x86), the sought value is in the bottom part */
156 if ( *(int8_t*)&checkval != 0 ) /*lint !e774*/
157 {
158 /* if the top part is nonzero, we assume that the number is negative */
159 if ( work.small[1] != 0 ) /*lint !e2662*/
160 {
161 work.big = -num;
162 return -work.small[0];
163 }
164 return work.small[0];
165 }
166
167 /* otherwise we have a big-endian machine (e.g., PowerPC); the sought value is in the top part */
168 assert( *(int8_t*)&checkval == 0 );
169
170 /* if the bottom part is nonzero, we assume that the number is negative */
171 if ( work.small[0] != 0 ) /*lint !e774*/
172 {
173 work.big = -num;
174 return -work.small[1]; /*lint !e2662*/
175 }
176 return work.small[1];
177}
178#endif
179
180/** returns whether Lapack s available, i.e., whether it has been linked in */
182 int* majorver, /**< major version number */
183 int* minorver, /**< minor version number */
184 int* patchver /**< patch version number */
185 )
186{
187#ifdef SCIP_WITH_LAPACK
188 LAPACKINTTYPE MAJOR = 0LL;
189 LAPACKINTTYPE MINOR = 0LL;
190 LAPACKINTTYPE PATCH = 0LL;
191#endif
192
193 assert( majorver != NULL );
194 assert( minorver != NULL );
195 assert( patchver != NULL );
196
197#ifdef SCIP_WITH_LAPACK
198 F77_FUNC(ilaver, ILAVER)(&MAJOR, &MINOR, &PATCH);
199
200 *majorver = convertToInt(MAJOR);
201 *minorver = convertToInt(MINOR);
202 *patchver = convertToInt(PATCH);
203#else
204 *majorver = -1;
205 *minorver = -1;
206 *patchver = -1;
207#endif
208}
209
210#ifdef SCIP_WITH_LAPACK
211/** computes eigenvalues of a symmetric matrix using LAPACK */
212static
213SCIP_RETCODE lapackComputeEigenvalues(
214 BMS_BUFMEM* bufmem, /**< buffer memory */
215 SCIP_Bool geteigenvectors, /**< Should also the eigenvectors be computed? */
216 int n, /**< size of matrix */
217 SCIP_Real* A, /**< matrix data on input (size n * n); eigenvectors on output if geteigenvectors == TRUE */
218 SCIP_Real* eigenvalues /**< pointer to store eigenvalue */
219 )
220{
221 LAPACKINTTYPE* IWORK;
222 LAPACKINTTYPE* ISUPPZ;
223 LAPACKINTTYPE N;
224 LAPACKINTTYPE INFO;
225 LAPACKINTTYPE LDA;
226 LAPACKINTTYPE WISIZE;
227 LAPACKINTTYPE IL;
228 LAPACKINTTYPE IU;
229 LAPACKINTTYPE M;
230 LAPACKINTTYPE LDZ;
231 LAPACKINTTYPE LWORK;
232 LAPACKINTTYPE LIWORK;
233 SCIP_Real* WORK;
234 SCIP_Real* WTMP;
235 SCIP_Real* Z = NULL;
236 SCIP_Real ABSTOL;
237 SCIP_Real WSIZE;
238 SCIP_Real VL;
239 SCIP_Real VU;
240 char JOBZ;
241 char RANGE;
242 char UPLO;
243
244 assert( bufmem != NULL );
245 assert( n > 0 );
246 assert( n < INT_MAX );
247 assert( A != NULL );
248 assert( eigenvalues != NULL );
249
250 N = n;
251 JOBZ = geteigenvectors ? 'V' : 'N';
252 RANGE = 'A';
253 UPLO = 'L';
254 LDA = n;
255 ABSTOL = 0.0; /* we use abstol = 0, since some lapack return an error otherwise */
256 VL = -1e20;
257 VU = 1e20;
258 IL = 0;
259 IU = n;
260 M = n;
261 LDZ = n;
262 INFO = 0LL;
263
264 /* standard LAPACK workspace query, to get the amount of needed memory */
265 LWORK = -1LL;
266 LIWORK = -1LL;
267
268 /* this computes the internally needed memory and returns this as (the first entries of [the 1x1 arrays]) WSIZE and WISIZE */
269 F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
270 &N, NULL, &LDA,
271 NULL, NULL,
272 &IL, &IU,
273 &ABSTOL, &M, NULL, NULL,
274 &LDZ, NULL, &WSIZE,
275 &LWORK, &WISIZE, &LIWORK,
276 &INFO);
277
278 if ( convertToInt(INFO) != 0 )
279 {
280 SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d.\n", convertToInt(INFO));
281 return SCIP_ERROR;
282 }
283
284 /* allocate workspace */
285 LWORK = SCIP_RealTOINT(WSIZE);
286 LIWORK = WISIZE;
287
288 SCIP_ALLOC( BMSallocBufferMemoryArray(bufmem, &WORK, (int) LWORK) );
289 SCIP_ALLOC( BMSallocBufferMemoryArray(bufmem, &IWORK, (int) LIWORK) );
290 SCIP_ALLOC( BMSallocBufferMemoryArray(bufmem, &WTMP, (int) N) );
291 SCIP_ALLOC( BMSallocBufferMemoryArray(bufmem, &ISUPPZ, 2) ); /*lint !e506*/
292 if ( geteigenvectors )
293 {
294 SCIP_ALLOC( BMSallocBufferMemoryArray(bufmem, &Z, n * n) ); /*lint !e647*/
295 }
296
297 /* call the function */
298 F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
299 &N, A, &LDA,
300 &VL, &VU,
301 &IL, &IU,
302 &ABSTOL, &M, WTMP, Z,
303 &LDZ, ISUPPZ, WORK,
304 &LWORK, IWORK, &LIWORK,
305 &INFO);
306
307 /* handle output */
308 if ( convertToInt(INFO) == 0 )
309 {
310 int m;
311 int i;
312 int j;
313
314 m = convertToInt(M);
315 for (i = 0; i < m; ++i)
316 eigenvalues[i] = WTMP[i];
317 for (i = m; i < n; ++i)
318 eigenvalues[i] = SCIP_INVALID;
319
320 /* possibly overwrite matrix with eigenvectors */
321 if ( geteigenvectors )
322 {
323 for (i = 0; i < m; ++i)
324 {
325 for (j = 0; j < n; ++j)
326 A[i * n + j] = Z[i * n + j];
327 }
328 }
329 }
330
331 /* free memory */
333 BMSfreeBufferMemoryArray(bufmem, &ISUPPZ);
334 BMSfreeBufferMemoryArray(bufmem, &WTMP);
335 BMSfreeBufferMemoryArray(bufmem, &IWORK);
336 BMSfreeBufferMemoryArray(bufmem, &WORK);
337
338 if ( convertToInt(INFO) != 0 )
339 {
340 SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d.\n", convertToInt(INFO));
341 return SCIP_ERROR;
342 }
343
344 return SCIP_OKAY;
345}
346#endif
347
348/** computes eigenvalues and eigenvectors of a dense symmetric matrix
349 *
350 * Calls Lapack's DSYEV function.
351 */
353 BMS_BUFMEM* bufmem, /**< buffer memory (or NULL if IPOPT is used) */
354 SCIP_Bool geteigenvectors, /**< should also eigenvectors should be computed? */
355 int N, /**< dimension */
356 SCIP_Real* a, /**< matrix data on input (size N*N); eigenvectors on output if geteigenvectors == TRUE */
357 SCIP_Real* w /**< array to store eigenvalues (size N) (or NULL) */
358 )
359{
360 /* if IPOPT is available, call its LAPACK routine */
362 {
363 SCIP_CALL( SCIPcallLapackDsyevIpopt(geteigenvectors, N, a, w) );
364 }
365 else
366 {
367 assert( bufmem != NULL );
368#ifdef SCIP_WITH_LAPACK
369 SCIP_CALL( lapackComputeEigenvalues(bufmem, geteigenvectors, N, a, w) );
370#else
371 SCIPerrorMessage("Lapack not available.\n");
372 return SCIP_PLUGINNOTFOUND;
373#endif
374 }
375
376 return SCIP_OKAY;
377}
378
379/** solves a linear problem of the form Ax = b for a regular matrix A
380 *
381 * Calls Lapacks DGETRF routine to calculate a LU factorization and uses this factorization to solve
382 * the linear problem Ax = b.
383 *
384 * Code taken from nlpi_ipopt.cpp
385 */
387 BMS_BUFMEM* bufmem, /**< buffer memory */
388 int n, /**< dimension */
389 SCIP_Real* A, /**< matrix data on input (size N*N); filled column-wise */
390 SCIP_Real* b, /**< right hand side vector (size N) */
391 SCIP_Real* x, /**< buffer to store solution (size N) */
392 SCIP_Bool* success /**< pointer to store if the solving routine was successful */
393 )
394{
395 assert( n > 0 );
396 assert( A != NULL );
397 assert( b != NULL );
398 assert( x != NULL );
399 assert( success != NULL );
400
401 /* if possible, use IPOPT */
403 {
404 SCIP_CALL( SCIPsolveLinearEquationsIpopt(n, A, b, x, success) );
405 }
406 else
407 {
408#ifdef SCIP_WITH_LAPACK
409 LAPACKINTTYPE INFO;
410 LAPACKINTTYPE N;
411 LAPACKINTTYPE* pivots;
412 SCIP_Real* Atmp = NULL;
413 SCIP_Real* btmp = NULL;
414
415 assert( bufmem != NULL );
416
417 SCIP_ALLOC( BMSduplicateBufferMemoryArray(bufmem, &Atmp, A, n * n) ); /*lint !e647*/
418 SCIP_ALLOC( BMSduplicateBufferMemoryArray(bufmem, &btmp, b, n) );
419 SCIP_ALLOC( BMSallocBufferMemoryArray(bufmem, &pivots, n) );
420
421 /* compute LU factorization */
422 N = n;
423 F77_FUNC(dgetrf, DGETRF)(&N, &N, Atmp, &N, pivots, &INFO);
424
425 if ( convertToInt(INFO) != 0 )
426 {
427 SCIPdebugMessage("There was an error when calling DGETRF. INFO = %d\n", convertToInt(INFO));
428 *success = FALSE;
429 }
430 else
431 {
432 LAPACKINTTYPE NRHS = 1LL;
433 char TRANS = 'N';
434
435 /* solve system */
436 F77_FUNC(dgetrs, DGETRS)(&TRANS, &N, &NRHS, Atmp, &N, pivots, btmp, &N, &INFO);
437
438 if ( convertToInt(INFO) != 0 )
439 {
440 SCIPdebugMessage("There was an error when calling DGETRF. INFO = %d\n", convertToInt(INFO));
441 *success = FALSE;
442 }
443 else
444 *success = TRUE;
445
446 /* copy the solution */
447 BMScopyMemoryArray(x, btmp, n);
448 }
449
450 BMSfreeBufferMemoryArray(bufmem, &pivots);
451 BMSfreeBufferMemoryArray(bufmem, &btmp);
452 BMSfreeBufferMemoryArray(bufmem, &Atmp);
453#else
454 SCIP_UNUSED(bufmem);
455
456 /* call fallback solution in nlpi_ipopt_dummy */
457 SCIP_CALL( SCIPsolveLinearEquationsIpopt(n, A, b, x, success) );
458#endif
459 }
460
461 return SCIP_OKAY;
462}
463
464/**@} */
SCIP_VAR * w
Definition: circlepacking.c:67
SCIP_VAR * a
Definition: circlepacking.c:66
SCIP_VAR ** b
Definition: circlepacking.c:65
SCIP_VAR ** x
Definition: circlepacking.c:63
common defines and data types used in all packages of SCIP
#define NULL
Definition: def.h:267
#define SCIP_UNUSED(x)
Definition: def.h:428
#define SCIP_INVALID
Definition: def.h:193
#define SCIP_Bool
Definition: def.h:91
#define SCIP_ALLOC(x)
Definition: def.h:385
#define SCIP_Real
Definition: def.h:173
#define TRUE
Definition: def.h:93
#define FALSE
Definition: def.h:94
#define SCIP_CALL(x)
Definition: def.h:374
SCIP_RETCODE SCIPcallLapackDsyevIpopt(SCIP_Bool computeeigenvectors, int N, SCIP_Real *a, SCIP_Real *w)
SCIP_RETCODE SCIPsolveLinearEquationsIpopt(int N, SCIP_Real *A, SCIP_Real *b, SCIP_Real *x, SCIP_Bool *success)
SCIP_Bool SCIPisIpoptAvailableIpopt(void)
SCIP_Bool SCIPlapackIsAvailable(void)
Definition: lapack_calls.c:121
SCIP_RETCODE SCIPlapackComputeEigenvalues(BMS_BUFMEM *bufmem, SCIP_Bool geteigenvectors, int N, SCIP_Real *a, SCIP_Real *w)
Definition: lapack_calls.c:352
SCIP_RETCODE SCIPlapackSolveLinearEquations(BMS_BUFMEM *bufmem, int n, SCIP_Real *A, SCIP_Real *b, SCIP_Real *x, SCIP_Bool *success)
Definition: lapack_calls.c:386
void SCIPlapackVersion(int *majorver, int *minorver, int *patchver)
Definition: lapack_calls.c:181
interface methods for lapack functions
memory allocation routines
#define BMSduplicateBufferMemoryArray(mem, ptr, source, num)
Definition: memory.h:737
#define BMSfreeBufferMemoryArray(mem, ptr)
Definition: memory.h:742
#define BMScopyMemoryArray(ptr, source, num)
Definition: memory.h:134
#define BMSallocBufferMemoryArray(mem, ptr, num)
Definition: memory.h:731
#define BMSfreeBufferMemoryArrayNull(mem, ptr)
Definition: memory.h:743
void F77_FUNC(filtersqp, FILTERSQP)
Ipopt NLP interface.
public methods for message output
#define SCIPerrorMessage
Definition: pub_message.h:64
#define SCIPdebugMessage
Definition: pub_message.h:96
type definitions for return codes for SCIP methods
@ SCIP_PLUGINNOTFOUND
Definition: type_retcode.h:54
@ SCIP_OKAY
Definition: type_retcode.h:42
@ SCIP_ERROR
Definition: type_retcode.h:43
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:63