SCIP-SDP  3.2.0
lapack_interface.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of SCIPSDP - a solving framework for mixed-integer */
4 /* semidefinite programs based on SCIP. */
5 /* */
6 /* Copyright (C) 2011-2013 Discrete Optimization, TU Darmstadt */
7 /* EDOM, FAU Erlangen-Nürnberg */
8 /* 2014-2020 Discrete Optimization, TU Darmstadt */
9 /* */
10 /* */
11 /* This program is free software; you can redistribute it and/or */
12 /* modify it under the terms of the GNU Lesser General Public License */
13 /* as published by the Free Software Foundation; either version 3 */
14 /* of the License, or (at your option) any later version. */
15 /* */
16 /* This program is distributed in the hope that it will be useful, */
17 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
18 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
19 /* GNU Lesser General Public License for more details. */
20 /* */
21 /* You should have received a copy of the GNU Lesser General Public License */
22 /* along with this program; if not, write to the Free Software */
23 /* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.*/
24 /* */
25 /* */
26 /* Based on SCIP - Solving Constraint Integer Programs */
27 /* Copyright (C) 2002-2020 Zuse Institute Berlin */
28 /* SCIP is distributed under the terms of the SCIP Academic Licence, */
29 /* see file COPYING in the SCIP distribution. */
30 /* */
31 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
32 
48 /* #define PRINTMATRICES /\* Should all matrices appearing in best rank-1 approximation heuristic be printed? *\/ */
49 
50 #include <assert.h>
51 
52 #include "lapack_interface.h"
53 #include "config.h" /* for F77_FUNC */
54 
55 #include "scip/def.h"
56 #include "scip/pub_message.h" /* for debug and error message */
57 #include "blockmemshell/memory.h"
58 #include "scip/type_retcode.h"
59 
60 /* turn off lint warnings for whole file: */
61 /*lint --e{788,818}*/
62 
63 
64 /* if we use 64 bit integers then use long long int, otherwise int */
65 #ifdef LAPACK_LONGLONGINT
66 typedef long long int LAPACKINTTYPE;
67 #else
68 typedef int LAPACKINTTYPE;
69 #endif
70 
71 
73 #define BMS_CALL(x) do \
74  { \
75  if( NULL == (x) ) \
76  { \
77  SCIPerrorMessage("No memory in function call\n"); \
78  return SCIP_NOMEMORY; \
79  } \
80  } \
81  while( FALSE )
82 
84 #define SCIP_RealTOINT(x) ((LAPACKINTTYPE) (x + 0.5))
85 
86 /*
87  * BLAS/LAPACK Calls
88  */
89 
94 void F77_FUNC(dsyevr, DSYEVR)(
95  char* JOBZ, char* RANGE, char* UPLO,
96  LAPACKINTTYPE* N, SCIP_Real* A, LAPACKINTTYPE* LDA,
97  SCIP_Real* VL, SCIP_Real* VU,
98  LAPACKINTTYPE* IL, LAPACKINTTYPE* IU,
99  SCIP_Real* ABSTOL, LAPACKINTTYPE* M, SCIP_Real* W, SCIP_Real* Z,
100  LAPACKINTTYPE* LDZ, LAPACKINTTYPE* ISUPPZ, SCIP_Real* WORK,
101  LAPACKINTTYPE* LWORK, LAPACKINTTYPE* IWORK, LAPACKINTTYPE* LIWORK,
102  int* INFO);
103 
104 
106 void F77_FUNC(dgemv, DGEMV)(char* TRANS, LAPACKINTTYPE* M,
107  LAPACKINTTYPE* N, SCIP_Real* ALPHA, SCIP_Real* A, LAPACKINTTYPE* LDA,
108  SCIP_Real* X, LAPACKINTTYPE* INCX, SCIP_Real* BETA, SCIP_Real* Y, LAPACKINTTYPE* INCY);
109 
110 
112 void F77_FUNC(dgemm, DGEMM)(char* TRANSA, char* TRANSB, LAPACKINTTYPE* M, LAPACKINTTYPE* N, LAPACKINTTYPE* K, SCIP_Real* ALPHA,
113  SCIP_Real* A, LAPACKINTTYPE* LDA, SCIP_Real* B, LAPACKINTTYPE* LDB, SCIP_Real* BETA, SCIP_Real* C, LAPACKINTTYPE* LDC );
114 
115 /* LAPACK Fortran subroutine DGELSD */
116 void F77_FUNC(dgelsd, DGELSD)(int* M, int* N, int* NRHS,
117  SCIP_Real* A, int* LDA, SCIP_Real* b, int* LDB, SCIP_Real* S, SCIP_Real* RCOND, int* RANK,
118  SCIP_Real* WORK, int* LWORK, int* IWORK, int* INFO );
119 
120 
124 /*
125  * Functions
126  */
127 
132 SCIP_RETCODE SCIPlapackComputeIthEigenvalue(
133  BMS_BUFMEM* bufmem,
134  SCIP_Bool geteigenvectors,
135  int n,
136  SCIP_Real* A,
137  int i,
138  SCIP_Real* eigenvalue,
139  SCIP_Real* eigenvector
140  )
141 {
142  LAPACKINTTYPE N;
143  LAPACKINTTYPE INFO;
144  char JOBZ;
145  char RANGE;
146  char UPLO;
147  LAPACKINTTYPE LDA;
148  SCIP_Real* WORK;
149  LAPACKINTTYPE LWORK;
150  LAPACKINTTYPE* IWORK;
151  LAPACKINTTYPE LIWORK;
152  SCIP_Real* WTMP;
153  SCIP_Real ABSTOL;
154  LAPACKINTTYPE IL;
155  LAPACKINTTYPE IU;
156  LAPACKINTTYPE M;
157  LAPACKINTTYPE LDZ;
158  SCIP_Real WSIZE;
159  LAPACKINTTYPE WISIZE;
160  SCIP_Real VL;
161  SCIP_Real VU;
162  LAPACKINTTYPE* ISUPPZ;
163 
164  assert( bufmem != NULL );
165  assert( n > 0 );
166  assert( A != NULL );
167  assert( 0 < i && i <= n );
168  assert( eigenvalue != NULL );
169  assert( ( ! geteigenvectors) || eigenvector != NULL );
170 
171  N = n;
172  JOBZ = geteigenvectors ? 'V' : 'N';
173  RANGE = 'I';
174  UPLO = 'L';
175  LDA = n;
176  ABSTOL = 0.0;
177  IL = i;
178  IU = i;
179  M = 1;
180  LDZ = n;
181 
182  /* standard LAPACK workspace query, to get the amount of needed memory */
183 #ifdef LAPACK_LONGLONGINT
184  LWORK = -1LL;
185  LIWORK = -1LL;
186 #else
187  LWORK = -1;
188  LIWORK = -1;
189 #endif
190 
191  /* this computes the internally needed memory and returns this as (the first entries of [the 1x1 arrays]) WSIZE and WISIZE */
192  F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
193  &N, NULL, &LDA,
194  NULL, NULL,
195  &IL, &IU,
196  &ABSTOL, &M, NULL, NULL,
197  &LDZ, NULL, &WSIZE,
198  &LWORK, &WISIZE, &LIWORK,
199  &INFO);
200 
201  /* for some reason this code seems to be called with INFO=0 within UG */
202  if ( INFO != 0 )
203  {
204  SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %lld.\n", INFO);
205  return SCIP_ERROR;
206  }
207 
208  /* allocate workspace */
209  LWORK = SCIP_RealTOINT(WSIZE);
210  LIWORK = WISIZE;
211 
212  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WORK, (int) LWORK) );
213  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &IWORK, (int) LIWORK) );
214  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WTMP, (int) N) );
215  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &ISUPPZ, 2) ); /*lint !e506*/
216 
217  /* call the function */
218  VL = -1e20;
219  VU = 1e20;
220 
221  F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
222  &N, A, &LDA,
223  &VL, &VU,
224  &IL, &IU,
225  &ABSTOL, &M, WTMP, eigenvector,
226  &LDZ, ISUPPZ, WORK,
227  &LWORK, IWORK, &LIWORK,
228  &INFO );
229 
230  if ( INFO != 0 )
231  {
232  SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d.\n", INFO);
233  return SCIP_ERROR;
234  }
235 
236  /* handle output */
237  *eigenvalue = WTMP[0];
238 
239  /* free memory */
240  BMSfreeBufferMemoryArray(bufmem, &ISUPPZ);
241  BMSfreeBufferMemoryArray(bufmem, &WTMP);/*lint !e737*/
242  BMSfreeBufferMemoryArray(bufmem, &IWORK);/*lint !e737*/
243  BMSfreeBufferMemoryArray(bufmem, &WORK);/*lint !e737*/
244 
245  return SCIP_OKAY;
246 }
247 
248 
250 SCIP_RETCODE SCIPlapackComputeEigenvectorDecomposition(
251  BMS_BUFMEM* bufmem,
252  int n,
253  SCIP_Real* A,
254  SCIP_Real* eigenvalues,
255  SCIP_Real* eigenvectors
256  )
257 {
258  LAPACKINTTYPE N;
259  LAPACKINTTYPE INFO;
260  char JOBZ;
261  char RANGE;
262  char UPLO;
263  LAPACKINTTYPE LDA;
264  SCIP_Real* WORK;
265  LAPACKINTTYPE LWORK;
266  LAPACKINTTYPE* IWORK;
267  LAPACKINTTYPE LIWORK;
268  SCIP_Real ABSTOL;
269  LAPACKINTTYPE M;
270  LAPACKINTTYPE LDZ;
271  SCIP_Real WSIZE;
272  LAPACKINTTYPE WISIZE;
273  SCIP_Real VL;
274  SCIP_Real VU;
275  LAPACKINTTYPE* ISUPPZ;
276 
277  assert( bufmem != NULL );
278  assert( n > 0 );
279  assert( A != NULL );
280  assert( eigenvalues != NULL );
281  assert( eigenvectors != NULL );
282 
283  N = n;
284  JOBZ = 'V';
285  RANGE = 'A';
286  UPLO = 'L';
287  LDA = n;
288  ABSTOL = 0.0;
289  M = n;
290  LDZ = n;
291 
292  /* standard LAPACK workspace query, to get the amount of needed memory */
293 #ifdef LAPACK_LONGLONGINT
294  LWORK = -1LL;
295  LIWORK = -1LL;
296 #else
297  LWORK = -1;
298  LIWORK = -1;
299 #endif
300 
301  /* this computes the internally needed memory and returns this as (the first entries of [the 1x1 arrays]) WSIZE and WISIZE */
302  F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
303  &N, NULL, &LDA,
304  NULL, NULL,
305  NULL, NULL,
306  &ABSTOL, &M, NULL, NULL,
307  &LDZ, NULL, &WSIZE,
308  &LWORK, &WISIZE, &LIWORK,
309  &INFO );
310 
311  if ( INFO != 0 )
312  {
313  SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d.\n", INFO);
314  return SCIP_ERROR;
315  }
316 
317  /* allocate workspace */
318  LWORK = SCIP_RealTOINT(WSIZE);
319  LIWORK = WISIZE;
320 
321  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WORK, (int) LWORK) );
322  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &IWORK, (int) LIWORK) );
323  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &ISUPPZ, (int) 2 * N) );
324 
325  /* call the function */
326  VL = -1e20;
327  VU = 1e20;
328 
329  F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
330  &N, A, &LDA,
331  &VL, &VU,
332  NULL, NULL,
333  &ABSTOL, &M, eigenvalues, eigenvectors,
334  &LDZ, ISUPPZ, WORK,
335  &LWORK, IWORK, &LIWORK,
336  &INFO );
337 
338  if ( INFO != 0 )
339  {
340  SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d.\n", INFO);
341  return SCIP_ERROR;
342  }
343 
344  /* free memory */
345  BMSfreeBufferMemoryArray(bufmem, &ISUPPZ);
346  BMSfreeBufferMemoryArray(bufmem, &IWORK);/*lint !e737*/
347  BMSfreeBufferMemoryArray(bufmem, &WORK);/*lint !e737*/
348 
349  return SCIP_OKAY;
350 }
351 
352 
354 SCIP_RETCODE SCIPlapackMatrixVectorMult(
355  int nrows,
356  int ncols,
357  SCIP_Real* matrix,
358  SCIP_Real* vector,
359  SCIP_Real* result
360  )
361 {
362  char TRANS;
363  LAPACKINTTYPE M;
364  LAPACKINTTYPE N;
365  SCIP_Real ALPHA;
366  SCIP_Real* A;
367  LAPACKINTTYPE LDA;
368  SCIP_Real* X;
369  LAPACKINTTYPE INCX;
370  SCIP_Real BETA;
371  SCIP_Real* Y;
372  LAPACKINTTYPE INCY;
373 
374  TRANS = 'N';
375  M = nrows;
376  N = ncols;
377  ALPHA = 1.0;
378  A = matrix;
379  LDA = nrows;
380  X = vector;
381  INCX = 1;
382  BETA = 0.0;
383  Y = result;
384  INCY = 1;
385 
386  F77_FUNC(dgemv, DGEMV)(&TRANS, &M, &N, &ALPHA, A, &LDA, X, &INCX, &BETA, Y, &INCY);
387 
388  return SCIP_OKAY;
389 }
390 
391 
393 SCIP_RETCODE SCIPlapackMatrixMatrixMult(
394  int nrowsA,
395  int ncolsA,
396  SCIP_Real* matrixA,
397  SCIP_Bool transposeA,
398  int nrowsB,
399  int ncolsB,
400  SCIP_Real* matrixB,
401  SCIP_Bool transposeB,
402  SCIP_Real* result
403  )
404 {
405  char TRANSA;
406  char TRANSB;
407  LAPACKINTTYPE M;
408  LAPACKINTTYPE N;
409  LAPACKINTTYPE K;
410  SCIP_Real ALPHA;
411  LAPACKINTTYPE LDA;
412  LAPACKINTTYPE LDB;
413  SCIP_Real BETA;
414  LAPACKINTTYPE LDC;
415 
416  assert( (transposeA && transposeB && (nrowsA == ncolsB)) || (transposeA && !transposeB && (nrowsA == nrowsB))
417  || (!transposeA && transposeB && (ncolsA == ncolsB)) || (!transposeA && !transposeB && (ncolsA == nrowsB)) );
418 
419  TRANSA = transposeA ? 'T' : 'N';
420  TRANSB = transposeB ? 'T' : 'N';
421  M = transposeA ? ncolsA : nrowsA;
422  N = transposeB ? nrowsB : ncolsB;
423  K = transposeA ? nrowsA : ncolsA;
424  ALPHA = 1.0;
425  LDA = transposeA ? K : M;
426  LDB = transposeB ? N : K;
427  BETA = 0.0;
428  LDC = M;
429 
430  F77_FUNC(dgemm, DGEMM)(&TRANSA, &TRANSB, &M, &N, &K, &ALPHA, matrixA, &LDA, matrixB, &LDB, &BETA, result, &LDC);
431 
432  return SCIP_OKAY;
433 }
434 
438 SCIP_RETCODE SCIPlapackLinearSolve(
439  BMS_BUFMEM* bufmem,
440  int m,
441  int n,
442  SCIP_Real* A,
443  SCIP_Real* b,
444  SCIP_Real* x
445  )
446 {
447  int i;
448  int M;
449  int N;
450  int NRHS;
451  int LDA;
452  int LDB;
453  SCIP_Real* S;
454  SCIP_Real RCOND;
455  int RANK;
456  SCIP_Real* WORK;
457  int LWORK;
458  int LIWORK;
459  int* IWORK;
460  int INFO;
461  SCIP_Real WSIZE;
462  int WISIZE;
463 #ifdef PRINTMATRICES
464  SCIP_Real residual = 0.0;
465 #endif
466 
467  assert( bufmem != NULL );
468  assert( m > 0 );
469  assert( n > 0 );
470  assert( A != NULL );
471  assert( b != NULL );
472  assert( x != NULL );
473 
474  M = m;
475  N = n;
476  NRHS = 1;
477  LDA = m;
478  LDB = MAX(M,N);
479  RCOND = 0.0;
480 
481  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &S, MIN(m,n)) );
482 
483  /* standard LAPACK workspace query, to get the amount of needed memory */
484  LWORK = -1LL;
485 
486  /* this computes the internally needed memory and returns this as (the first entry of [the 1x1 array]) WSIZE */
487  F77_FUNC(dgelsd, DGELSD)( &M, &N, &NRHS,
488  NULL, &LDA, NULL, &LDB, NULL,
489  &RCOND, &RANK, &WSIZE, &LWORK,
490  &WISIZE, &INFO );
491 
492  if ( INFO != 0 )
493  {
494  SCIPerrorMessage("There was an error when calling DGELSD. INFO = %d\n", INFO);
495  return SCIP_ERROR;
496  }
497 
498  /* allocate workspace */
499  LWORK = SCIP_RealTOINT(WSIZE);
500  LIWORK = WISIZE;
501 
502  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WORK, (int) LWORK) );
503  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &IWORK, (int) LIWORK) );
504 
505  /* call the function */
506  F77_FUNC(dgelsd, DGELSD)( &M, &N, &NRHS,
507  A, &LDA, b, &LDB, S, &RCOND, &RANK,
508  WORK, &LWORK, IWORK, &INFO );
509 
510 #ifdef PRINTMATRICES
511  printf("LWORK = %d\n", LWORK);
512  printf("LIWORK = %d\n", LIWORK);
513  printf("A has size (%d,%d), is of rank %d\n", M, N, RANK);
514  printf("Minimum l2-norm solution of linear equation system:\n");
515 
516  for (i = 0; i < n; ++i)
517  {
518  printf("(%d, %f) ", i, b[i]);
519  }
520 
521  for (i = n; i < m; ++i)
522  {
523  residual += b[i] * b[i];
524  }
525  printf("\n");
526 
527  printf("Residual sum-of-squares for the solution is %f\n", residual);
528 #endif
529 
530  if ( INFO != 0 )
531  {
532  SCIPerrorMessage("There was an error when calling DGELSD. INFO = %d\n", INFO);
533  return SCIP_ERROR;
534  }
535 
536  /* LAPACK overwrites the right-hand side with the result */
537  for (i = 0; i < n; ++i)
538  x[i] = b[i];
539 
540  /* free memory */
541  BMSfreeBufferMemoryArray(bufmem, &IWORK);/*lint !e737*/
542  BMSfreeBufferMemoryArray(bufmem, &WORK);/*lint !e737*/
543  BMSfreeBufferMemoryArray(bufmem, &S);/*lint !e737*/
544 
545  return SCIP_OKAY;
546 }
547 
SCIP_RealTOINT
#define SCIP_RealTOINT(x)
Definition: lapack_interface.c:84
config.h
Macros for calling Fortran-functions.
BMS_CALL
#define BMS_CALL(x)
Definition: lapack_interface.c:73
SCIPlapackMatrixMatrixMult
SCIP_RETCODE SCIPlapackMatrixMatrixMult(int nrowsA, int ncolsA, SCIP_Real *matrixA, SCIP_Bool transposeA, int nrowsB, int ncolsB, SCIP_Real *matrixB, SCIP_Bool transposeB, SCIP_Real *result)
Definition: lapack_interface.c:393
SCIPlapackMatrixVectorMult
SCIP_RETCODE SCIPlapackMatrixVectorMult(int nrows, int ncols, SCIP_Real *matrix, SCIP_Real *vector, SCIP_Real *result)
Definition: lapack_interface.c:354
SCIPlapackComputeIthEigenvalue
SCIP_EXPORT SCIP_RETCODE SCIPlapackComputeIthEigenvalue(BMS_BUFMEM *bufmem, SCIP_Bool geteigenvectors, int n, SCIP_Real *A, int i, SCIP_Real *eigenvalue, SCIP_Real *eigenvector)
LAPACKINTTYPE
int LAPACKINTTYPE
Definition: lapack_interface.c:68
SCIPlapackLinearSolve
SCIP_RETCODE SCIPlapackLinearSolve(BMS_BUFMEM *bufmem, int m, int n, SCIP_Real *A, SCIP_Real *b, SCIP_Real *x)
Definition: lapack_interface.c:438
lapack_interface.h
interface methods for eigenvector computation and matrix multiplication using openblas ...
F77_FUNC
#define F77_FUNC(name, NAME)
Definition: config.h:43
SCIPlapackComputeEigenvectorDecomposition
SCIP_RETCODE SCIPlapackComputeEigenvectorDecomposition(BMS_BUFMEM *bufmem, int n, SCIP_Real *A, SCIP_Real *eigenvalues, SCIP_Real *eigenvectors)
Definition: lapack_interface.c:250
Generated by   doxygen 1.8.11