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