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