source: trunk/Couenne/src/cut/sdpcuts/dsyevx_wrapper.cpp @ 944

Last change on this file since 944 was 944, checked in by pbelotti, 8 years ago

switched from dsyevx to dsyev as that is used in coinlapack

  • Property svn:keywords set to Author Date Id Revision
File size: 2.6 KB
Line 
1/* $Id: dsyevx_wrapper.cpp 944 2013-03-29 03:36:09Z pbelotti $
2 *
3 * Name:    dsyevx_rapper.cpp
4 * Authors: Andrea Qualizza
5 *          Pietro Belotti
6 * Purpose:
7 *
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11#include <stdlib.h>
12#include <stdio.h>
13#include <math.h>
14
15#include "CoinFinite.hpp"
16#include "CouenneConfig.h"
17
18//#define DEBUG
19
20extern "C" {
21
22  /* Lapack routine to compute orthonormal eigenvalues/eigenvectors (in Fortran) */
23
24  void F77_FUNC(dsyev,DSYEV) (
25                                  char   *,
26                                  char   *,
27                                  char   *,
28                                  int    *,
29                                  double *,
30                                  int    *,
31                                  double *,
32                                  double *,
33                                  int    *,
34                                  int    *,
35                                  double *,
36                                  int    *,
37                                  double *,
38                                  double *,
39                                  int    *,
40                                  double *,
41                                  int    *,
42                                  int    *,
43                                  int    *,
44                                  int    *);
45}
46
47
48int dsyevx_interface (int n, double *A, int &m, 
49                      double * &w, 
50                      double * &z, // output values
51                      double tolerance,
52                      double lb_ev, 
53                      double ub_ev,
54                      int firstidx,
55                      int lastidx) {
56
57#ifdef DEBUG
58
59  printf ("matrix:\n---------------------------------\n");
60  for (int   i=0; i<n; ++i) {
61    for (int j=0; j<n; ++j)
62      printf ("%g ", A [i*n+j]);
63    printf ("\n");
64  }
65  printf ("---------------------------------\n");
66#endif
67
68  if (NULL == w) w = new double [n];
69  if (NULL == z) z = new double [n*n];
70
71  m = n;
72
73  int lwork = 8*n;
74
75  char jobz  = 'V';  // compute both eigenvalues and eigenvectors
76  char range = 'V';  // range for selection is on values of eigenvalues
77  char uplo  = 'U';  // upper triangular matrix is given
78
79  int il  = firstidx; // index of the eigenvalue to be returned 1=first
80  int iu  = lastidx;  // index of the last eigenvalue to be returuned 1=first
81  int lda = n;        // leading dimension of A
82  int ldz = n;        // leading dimension of z
83
84  int info; // output status
85
86  int *ifail = new int [n];
87  int *iwork = new int [5*n]; 
88 
89  double abstol = tolerance;    // absolute tolerance
90  double vl     = lb_ev;        // minimum eigenvalue wanted
91  double vu     = ub_ev;        // maximum
92
93  double *work  = new double [lwork];
94
95  // Equivalent:
96  // Ipopt::IpLapackDsyev (true, n, A, lda, w, info);
97
98  F77_FUNC
99    (dsyev,DSYEV)
100    (&jobz, &range, &uplo, &n, 
101     A, &lda, 
102     &vl, &vu, &il, &iu,
103     &abstol, &m, 
104     w, z, &ldz, work, &lwork, iwork, ifail, &info);
105
106  if (info) {
107    printf (":: dsyevx returned status %d\n", info);
108#ifdef CHECK
109    for(int i=0; i<m; i++) {
110      if(ifail[i] > 0) {
111        printf("### WARNING: dsyevx_wrapper(): ifail[%d]: %d   curr_ev[%d]=%.18f\n"
112               , i, ifail [i], ifail [i], w [ifail [i]]);
113      }
114    }
115#endif
116  }
117
118  delete [] work;
119  delete [] ifail;
120  delete [] iwork;
121
122  return m;
123}
Note: See TracBrowser for help on using the repository browser.