1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the program and library */ 4 /* SCIP --- Solving Constraint Integer Programs */ 5 /* */ 6 /* Copyright 2002-2022 Zuse Institute Berlin */ 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 */ 45 #include "blockmemshell/memory.h" 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 */ 55 typedef 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 */ 89 void 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 */ 99 void F77_FUNC(dgetrf, DGETRF)(LAPACKINTTYPE* M, LAPACKINTTYPE* N, SCIP_Real* A, 100 LAPACKINTTYPE* LDA, LAPACKINTTYPE* IPIV, LAPACKINTTYPE* INFO); 101 102 /** LAPACK Fortran subroutine DGETRS */ 103 void 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 */ 108 void 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 */ 121 SCIP_Bool SCIPlapackIsAvailable(void) 122 { 123 if ( SCIPisIpoptAvailableIpopt() ) 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 */ 139 static 140 int 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 */ 181 void SCIPlapackVersion( 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 */ 212 static 213 SCIP_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 */ 332 BMSfreeBufferMemoryArrayNull(bufmem, &Z); 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 */ 352 SCIP_RETCODE SCIPlapackComputeEigenvalues( 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 */ 361 if ( SCIPisIpoptAvailableIpopt() ) 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 */ 386 SCIP_RETCODE SCIPlapackSolveLinearEquations( 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 */ 402 if ( SCIPisIpoptAvailableIpopt() ) 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 /**@} */ 465