SCIP-SDP  3.1.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
lapack_sdpa.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-2018 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-2018 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 
33 /* #define SCIP_DEBUG*/
34 /* #define SCIP_MORE_DEBUG*/
35 
43 #include <assert.h>
44 
45 #include "lapack.h"
46 #include "config.h" /* for F77_FUNC */
47 
48 #include "scip/def.h"
49 #include "scip/pub_message.h" /* for debug and error message */
50 #include "blockmemshell/memory.h"
51 #include "scip/type_retcode.h"
52 
53 /* turn off lint warnings for whole file: */
54 /*lint --e{788,818}*/
55 
56 typedef long long int LAPACKINTTYPE;
57 #define SDPA_VERSION 740
58 
60 #define BMS_CALL(x) do \
61  { \
62  if( NULL == (x) ) \
63  { \
64  SCIPerrorMessage("No memory in function call\n"); \
65  return SCIP_NOMEMORY; \
66  } \
67  } \
68  while( FALSE )
69 
71 #define SCIP_RealTOINT(x) ((LAPACKINTTYPE) (x + 0.5))
72 
73 /*
74  * BLAS/LAPACK Calls
75  */
76 
81 void F77_FUNC(dsyevr, DSYEVR)(
82  char* JOBZ, char* RANGE, char* UPLO,
83  LAPACKINTTYPE* N, SCIP_Real* A, LAPACKINTTYPE* LDA,
84  SCIP_Real* VL, SCIP_Real* VU,
85  LAPACKINTTYPE* IL, LAPACKINTTYPE* IU,
86  SCIP_Real* ABSTOL, LAPACKINTTYPE* M, SCIP_Real* W, SCIP_Real* Z,
87  LAPACKINTTYPE* LDZ, LAPACKINTTYPE* ISUPPZ, SCIP_Real* WORK,
88  LAPACKINTTYPE* LWORK, LAPACKINTTYPE* IWORK, LAPACKINTTYPE* LIWORK,
89  int* INFO );
90 
91 
93 void F77_FUNC(dgemv, DGEMV)(char* TRANS, LAPACKINTTYPE* M,
94  LAPACKINTTYPE* N, SCIP_Real* ALPHA, SCIP_Real* A, LAPACKINTTYPE* LDA,
95  SCIP_Real* X, LAPACKINTTYPE* INCX, SCIP_Real* BETA, SCIP_Real* Y, LAPACKINTTYPE* INCY);
96 
97 
99 void F77_FUNC(dgemm, DGEMM)(char* TRANSA, char* TRANSB, LAPACKINTTYPE* M, LAPACKINTTYPE* N, LAPACKINTTYPE* K, SCIP_Real* ALPHA,
100  SCIP_Real* A, LAPACKINTTYPE* LDA, SCIP_Real* B, LAPACKINTTYPE* LDB, SCIP_Real* BETA, SCIP_Real* C, LAPACKINTTYPE* LDC );
101 
102 
106 /*
107  * Functions
108  */
109 
114 SCIP_RETCODE SCIPlapackComputeIthEigenvalue(
115  BMS_BUFMEM* bufmem,
116  SCIP_Bool geteigenvectors,
117  int n,
118  SCIP_Real* A,
119  int i,
120  SCIP_Real* eigenvalue,
121  SCIP_Real* eigenvector
122  )
123 {
124  LAPACKINTTYPE N;
125 #if ( SDPA_VERSION == 740 )
126  LAPACKINTTYPE INFO;
127 #else
128  int INFO;
129 #endif
130  char JOBZ;
131  char RANGE;
132  char UPLO;
133  LAPACKINTTYPE LDA;
134  SCIP_Real* WORK;
135  LAPACKINTTYPE LWORK;
136  LAPACKINTTYPE* IWORK;
137  LAPACKINTTYPE LIWORK;
138  SCIP_Real* WTMP;
139  SCIP_Real ABSTOL;
140  LAPACKINTTYPE IL;
141  LAPACKINTTYPE IU;
142  LAPACKINTTYPE M;
143  LAPACKINTTYPE LDZ;
144  SCIP_Real WSIZE;
145  LAPACKINTTYPE WISIZE;
146  SCIP_Real VL;
147  SCIP_Real VU;
148  LAPACKINTTYPE* ISUPPZ;
149 
150  assert( bufmem != NULL );
151  assert( n > 0 );
152  assert( A != NULL );
153  assert( 0 < i && i <= n );
154  assert( eigenvalue != NULL );
155  assert( ( ! geteigenvectors) || eigenvector != NULL );
156 
157  N = n;
158  JOBZ = geteigenvectors ? 'V' : 'N';
159  RANGE = 'I';
160  UPLO = 'L';
161  LDA = n;
162  ABSTOL = 0.0;
163  IL = i;
164  IU = i;
165  M = 1;
166  LDZ = n;
167 
168  /* standard LAPACK workspace query, to get the amount of needed memory */
169  LWORK = -1LL;
170  LIWORK = -1LL;
171 
172  /* this computes the internally needed memory and returns this as (the first entries of [the 1x1 arrays]) WSIZE and WISIZE */
173  F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
174  &N, NULL, &LDA,
175  NULL, NULL,
176  &IL, &IU,
177  &ABSTOL, &M, NULL, NULL,
178  &LDZ, NULL, &WSIZE,
179  &LWORK, &WISIZE, &LIWORK,
180  &INFO );
181 /* for some reason this code seems to be called with INFO=0 within UG */
182  if ( INFO != 0 )
183  {
184  SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %lld\n", INFO);
185  return SCIP_ERROR;
186  }
187 
188  /* allocate workspace */
189  LWORK = SCIP_RealTOINT(WSIZE);
190  LIWORK = WISIZE;
191 
192  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WORK, (int) LWORK) );
193  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &IWORK, (int) LIWORK) );
194  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WTMP, (int) N) );
195  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &ISUPPZ, 2) ); /*lint !e506*/
196 
197  /* call the function */
198  VL = -1e20;
199  VU = 1e20;
200 
201  F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
202  &N, A, &LDA,
203  &VL, &VU,
204  &IL, &IU,
205  &ABSTOL, &M, WTMP, eigenvector,
206  &LDZ, ISUPPZ, WORK,
207  &LWORK, IWORK, &LIWORK,
208  &INFO );
209 
210  if ( INFO != 0 )
211  {
212  SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d\n", INFO);
213  return SCIP_ERROR;
214  }
215 
216  /* handle output */
217  *eigenvalue = WTMP[0];
218 
219  /* free memory */
220  BMSfreeBufferMemoryArray(bufmem, &ISUPPZ);
221  BMSfreeBufferMemoryArray(bufmem, &WTMP);/*lint !e737*/
222  BMSfreeBufferMemoryArray(bufmem, &IWORK);/*lint !e737*/
223  BMSfreeBufferMemoryArray(bufmem, &WORK);/*lint !e737*/
224 
225  return SCIP_OKAY;
226 }
227 
228 
230 SCIP_RETCODE SCIPlapackComputeEigenvectorDecomposition(
231  BMS_BUFMEM* bufmem,
232  int n,
233  SCIP_Real* A,
234  SCIP_Real* eigenvalues,
235  SCIP_Real* eigenvectors
236  )
237 {
238  LAPACKINTTYPE N;
239 #if ( SDPA_VERSION == 740 )
240  LAPACKINTTYPE INFO;
241 #else
242  int INFO;
243 #endif
244  char JOBZ;
245  char RANGE;
246  char UPLO;
247  LAPACKINTTYPE LDA;
248  SCIP_Real* WORK;
249  LAPACKINTTYPE LWORK;
250  LAPACKINTTYPE* IWORK;
251  LAPACKINTTYPE LIWORK;
252  SCIP_Real ABSTOL;
253  LAPACKINTTYPE M;
254  LAPACKINTTYPE LDZ;
255  SCIP_Real WSIZE;
256  LAPACKINTTYPE WISIZE;
257  SCIP_Real VL;
258  SCIP_Real VU;
259  LAPACKINTTYPE* ISUPPZ;
260 
261  assert( bufmem != NULL );
262  assert( n > 0 );
263  assert( A != NULL );
264  assert( eigenvalues != NULL );
265  assert( eigenvectors != NULL );
266 
267  N = n;
268  JOBZ = 'V';
269  RANGE = 'A';
270  UPLO = 'L';
271  LDA = n;
272  ABSTOL = 0.0;
273  M = n;
274  LDZ = n;
275 
276  /* standard LAPACK workspace query, to get the amount of needed memory */
277  LWORK = -1LL;
278  LIWORK = -1LL;
279 
280  /* this computes the internally needed memory and returns this as (the first entries of [the 1x1 arrays]) WSIZE and WISIZE */
281  F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
282  &N, NULL, &LDA,
283  NULL, NULL,
284  NULL, NULL,
285  &ABSTOL, &M, NULL, NULL,
286  &LDZ, NULL, &WSIZE,
287  &LWORK, &WISIZE, &LIWORK,
288  &INFO );
289 
290  if ( INFO != 0 )
291  {
292  SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d\n", INFO);
293  return SCIP_ERROR;
294  }
295 
296  /* allocate workspace */
297  LWORK = SCIP_RealTOINT(WSIZE);
298  LIWORK = WISIZE;
299 
300  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WORK, (int) LWORK) );
301  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &IWORK, (int) LIWORK) );
302  BMS_CALL( BMSallocBufferMemoryArray(bufmem, &ISUPPZ, (int) 2 * N) );
303 
304  /* call the function */
305  VL = -1e20;
306  VU = 1e20;
307 
308  F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
309  &N, A, &LDA,
310  &VL, &VU,
311  NULL, NULL,
312  &ABSTOL, &M, eigenvalues, eigenvectors,
313  &LDZ, ISUPPZ, WORK,
314  &LWORK, IWORK, &LIWORK,
315  &INFO );
316 
317  if ( INFO != 0 )
318  {
319  SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d\n", INFO);
320  return SCIP_ERROR;
321  }
322 
323  /* free memory */
324  BMSfreeBufferMemoryArray(bufmem, &ISUPPZ);
325  BMSfreeBufferMemoryArray(bufmem, &IWORK);/*lint !e737*/
326  BMSfreeBufferMemoryArray(bufmem, &WORK);/*lint !e737*/
327 
328  return SCIP_OKAY;
329 }
330 
331 
333 SCIP_RETCODE SCIPlapackMatrixVectorMult(
334  int nrows,
335  int ncols,
336  SCIP_Real* matrix,
337  SCIP_Real* vector,
338  SCIP_Real* result
339  )
340 {
341  char TRANS;
342  LAPACKINTTYPE M;
343  LAPACKINTTYPE N;
344  SCIP_Real ALPHA;
345  SCIP_Real* A;
346  LAPACKINTTYPE LDA;
347  SCIP_Real* X;
348  LAPACKINTTYPE INCX;
349  SCIP_Real BETA;
350  SCIP_Real* Y;
351  LAPACKINTTYPE INCY;
352 
353  TRANS = 'N';
354  M = nrows;
355  N = ncols;
356  ALPHA = 1.0;
357  A = matrix;
358  LDA = nrows;
359  X = vector;
360  INCX = 1;
361  BETA = 0.0;
362  Y = result;
363  INCY = 1;
364 
365  F77_FUNC(dgemv, DGEMV)(&TRANS, &M, &N, &ALPHA, A, &LDA, X, &INCX, &BETA, Y, &INCY);
366 
367  return SCIP_OKAY;
368 }
369 
371 SCIP_RETCODE SCIPlapackMatrixMatrixMult(
372  int nrowsA,
373  int ncolsA,
374  SCIP_Real* matrixA,
375  SCIP_Bool transposeA,
376  int nrowsB,
377  int ncolsB,
378  SCIP_Real* matrixB,
379  SCIP_Bool transposeB,
380  SCIP_Real* result
381  )
382 {
383  char TRANSA;
384  char TRANSB;
385  LAPACKINTTYPE M;
386  LAPACKINTTYPE N;
387  LAPACKINTTYPE K;
388  SCIP_Real ALPHA;
389  LAPACKINTTYPE LDA;
390  LAPACKINTTYPE LDB;
391  SCIP_Real BETA;
392  LAPACKINTTYPE LDC;
393 
394  assert( (transposeA && transposeB && (nrowsA == ncolsB)) || (transposeA && !transposeB && (nrowsA == nrowsB))
395  || (!transposeA && transposeB && (ncolsA == ncolsB)) || (!transposeA && !transposeB && (ncolsA == nrowsB)) );
396 
397  TRANSA = transposeA ? 'T' : 'N';
398  TRANSB = transposeB ? 'T' : 'N';
399  M = transposeA ? ncolsA : nrowsA;
400  N = transposeB ? nrowsB : ncolsB;
401  K = transposeA ? nrowsA : ncolsA;
402  ALPHA = 1.0;
403  LDA = transposeA ? K : M;
404  LDB = transposeB ? N : K;
405  BETA = 0.0;
406  LDC = M;
407 
408  F77_FUNC(dgemm, DGEMM)(&TRANSA, &TRANSB, &M, &N, &K, &ALPHA, matrixA, &LDA, matrixB, &LDB, &BETA, result, &LDC);
409 
410  return SCIP_OKAY;
411 }
412 
LAPACKINTTYPE
long long int LAPACKINTTYPE
Definition: lapack_sdpa.c:56
SCIPlapackComputeEigenvectorDecomposition
SCIP_RETCODE SCIPlapackComputeEigenvectorDecomposition(BMS_BUFMEM *bufmem, int n, SCIP_Real *A, SCIP_Real *eigenvalues, SCIP_Real *eigenvectors)
Definition: lapack_sdpa.c:230
SCIPlapackComputeIthEigenvalue
EXTERN SCIP_RETCODE SCIPlapackComputeIthEigenvalue(BMS_BUFMEM *bufmem, SCIP_Bool geteigenvectors, int n, SCIP_Real *A, int i, SCIP_Real *eigenvalue, SCIP_Real *eigenvector)
config.h
Macros for calling Fortran-functions.
SCIPlapackMatrixVectorMult
SCIP_RETCODE SCIPlapackMatrixVectorMult(int nrows, int ncols, SCIP_Real *matrix, SCIP_Real *vector, SCIP_Real *result)
Definition: lapack_sdpa.c:333
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_sdpa.c:371
BMS_CALL
#define BMS_CALL(x)
Definition: lapack_sdpa.c:60
SCIP_RealTOINT
#define SCIP_RealTOINT(x)
Definition: lapack_sdpa.c:71
F77_FUNC
#define F77_FUNC(name, NAME)
Definition: config.h:43
lapack.h
interface methods for eigenvector computation and matrix multiplication using different versions of L...
Generated on Fri Jun 29 2018 06:53:24 for SCIP-SDP by   doxygen 1.8.7