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