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