/home/coin/SVN-release/OS-2.4.2/Couenne/src/cut/sdpcuts/dsyevx_wrapper.cpp

Go to the documentation of this file.
00001 /* $Id: dsyevx_wrapper.cpp 508 2011-02-15 21:52:44Z pbelotti $
00002  *
00003  * Name:    dsyevx_rapper.cpp
00004  * Author:  Andrea Qualizza
00005  * Purpose: 
00006  *
00007  * This file is licensed under the Eclipse Public License (EPL)
00008  */
00009 
00010 #include <stdlib.h>
00011 #include <stdio.h>
00012 #include <math.h>
00013 
00014 #include <misc_util.hpp>
00015 
00016 #include "CoinHelperFunctions.hpp"
00017 
00018 #include <tracer.hpp>
00019 
00020 #define ABS_TOL_EIG 1e-15
00021 
00022 
00023 
00024 void _dsyevx_value_range_wrapper (int n, double *A, int &m, double * &w, double * &z,double tolerance,double lb_ev, double ub_ev);
00025 
00026 void _dsyevx_index_range_wrapper (int n, double *A, int &m, double * &w, double * &z, double tolerance, int firstidx,int lastidx);
00027 
00028 
00029 #if 0
00030 extern "C" {
00031 
00032 /* Lapack routine to compute orthonormal eigenvalues/eigenvectors (in Fortran) */
00033 void dsyevx_ (char   *,
00034               char   *,
00035               char   *,
00036               int    *,
00037               double *,
00038               int    *,
00039               double *,
00040               double *,
00041               int    *,
00042               int    *,
00043               double *,
00044               int    *,
00045               double *,
00046               double *,
00047               int    *,
00048               double *,
00049               int    *,
00050               int    *,
00051               int    *,
00052               int    *);
00053 }
00054 
00055 #else
00056 
00057 /* Lapack routine to compute orthonormal eigenvalues/eigenvectors (in Fortran) */
00058 void dsyevx_ (char   *,
00059               char   *,
00060               char   *,
00061               int    *,
00062               double *,
00063               int    *,
00064               double *,
00065               double *,
00066               int    *,
00067               int    *,
00068               double *,
00069               int    *,
00070               double *,
00071               double *,
00072               int    *,
00073               double *,
00074               int    *,
00075               int    *,
00076               int    *,
00077               int    *) {}
00078 #endif
00079 
00080 
00081 
00082 void dsyevx_full_wrapper (int n, double *A, int &m, double * &w, double * &z, Tracer *tracer) {
00083         tracer->incrementMainTotalEigendecompositions();
00084         _dsyevx_value_range_wrapper (n,A,m,w,z,ABS_TOL_EIG,-COIN_DBL_MAX,COIN_DBL_MAX);
00085 }
00086 
00087 void dsyevx_wrapper_only_positive (int n, double *A, int &m, double * &w, double * &z, Tracer *tracer) {
00088         tracer->incrementMainTotalEigendecompositions();
00089         _dsyevx_value_range_wrapper (n,A,m,w,z,ABS_TOL_EIG,0.0,COIN_DBL_MAX);
00090 }
00091 
00092 void dsyevx_wrapper_only_negative (int n, double *A, int &m, double * &w, double * &z, Tracer *tracer) {
00093         tracer->incrementMainTotalEigendecompositions();
00094         _dsyevx_value_range_wrapper (n,A,m,w,z,ABS_TOL_EIG,-COIN_DBL_MAX,ABS_TOL_EIG);
00095 }
00096 
00097 void dsyevx_wrapper_only_most_neg (int n, double *A, int &m, double * &w, double * &z, Tracer *tracer) {
00098         tracer->incrementMainTotalEigendecompositions();
00099         _dsyevx_index_range_wrapper (n,A,m,w,z,ABS_TOL_EIG,1,1);
00100 }
00101 
00102 void dsyevx_wrapper_first_p (int n, double *A, int &m, double * &w, double * &z, int p, Tracer *tracer) {
00103         tracer->incrementMainTotalEigendecompositions();
00104         _dsyevx_index_range_wrapper (n,A,m,w,z,ABS_TOL_EIG,1,p);
00105 }
00106 
00107 //########################################################################################
00108 
00109 
00110 void _dsyevx_value_range_wrapper (int n, double *A, int &m, double * &w, double * &z,double tolerance,double lb_ev, double ub_ev) {
00111 // the lapack call destroys the original matrix A
00112         if(w == NULL)
00113                 w = new double[n];
00114         if(z == NULL)
00115                 z = new double[n*n];
00116 
00117         m = n;
00118 
00119         int lwork = 8*n;
00120 
00121         char jobz  = 'V';  // compute both eigenvalues and eigenvectors
00122         char range = 'V';  // range for selection is on values of eigenvalues
00123         char uplo  = 'U';  // upper triangular matrix is given
00124 
00125         int il     = 1;   // first  eigenvalue to be returned (not used)
00126         int iu     = n;   // second                           (not used)
00127         int info;         // output status
00128         int lda    = n;   // leading dimension of A
00129         int ldz    = n;   // leading dimension of z
00130 
00131         int *ifail = new int[n];
00132         int *iwork = new int[5*n]; 
00133 
00134  
00135         double abstol = tolerance;      // absolute tolerance
00136         double vl     = lb_ev;          // minimum eigenvalue wanted
00137         double vu     = ub_ev;          // maximum
00138 
00139         double *work  = new double[lwork];
00140 
00141         dsyevx_ (&jobz, &range, &uplo, &n, A, &lda, &vl, &vu, &il, &iu,
00142                 &abstol, &m, w, z, &ldz, work, &lwork, iwork, ifail, &info);
00143 
00144         if (info) {
00145                 printf (":: dsyevx returned status %d\n", info);
00146 #ifdef CHECK
00147                 for(int i=0; i<m; i++) {
00148                         if(ifail[i] > 0) {
00149                                 printf("### WARNING: dsyevx_wrapper(): ifail[%d]: %d   curr_ev[%d]=%.18f\n"
00150                                         ,i, ifail[i],ifail[i],w[ifail[i]]);
00151                         }
00152                 }
00153 #endif
00154         }
00155 
00156         if ((info == 0) && (m > 0)) { // there is at least one eigenvector
00157 #ifdef TRACE_ALL
00158                 cpp_printvecDBL("eigenvalues", w, m);
00159 #endif
00160         }
00161         
00162         delete [] work;
00163         delete [] ifail;
00164         delete [] iwork;
00165 }
00166 
00167 
00168 //########################################################################################
00169 
00170 
00171 void _dsyevx_index_range_wrapper (int n, double *A, int &m, double * &w, double * &z, double tolerance, int firstidx,int lastidx) {
00172 // the lapack call destroys the original matrix A
00173         if(w == NULL)
00174                 w = new double[n];
00175         if(z == NULL)
00176                 z = new double[n*n];
00177 
00178         m = n;
00179 
00180         int lwork = 8*n;
00181 
00182         char jobz  = 'V';  // compute both eigenvalues and eigenvectors
00183         char range = 'I';  // range selection
00184         char uplo  = 'U';  // upper triangular matrix is given
00185 
00186         int il     = firstidx;   // index of the eigenvalue to be returned 1=first
00187         int iu     = lastidx;   // index of the last eigenvalue to be returuned 1=first
00188         int info;         // output status
00189         int lda    = n;   // leading dimension of A
00190         int ldz    = n;   // leading dimension of z
00191 
00192         int *ifail = new int[n];
00193         int *iwork = new int[5*n]; 
00194 
00195  
00196         double abstol = ABS_TOL_EIG;    // absolute tolerance
00197         double vl     = 0.0;            // not used
00198         double vu     = 0.0;            // not used
00199 
00200         double *work  = new double[lwork];
00201 
00202         dsyevx_ (&jobz, &range, &uplo, &n, A, &lda, &vl, &vu, &il, &iu,
00203                 &abstol, &m, w, z, &ldz, work, &lwork, iwork, ifail, &info);
00204 
00205         if (info) {
00206                 printf (":: dsyevx returned status %d\n", info);
00207 #ifdef CHECK
00208                 for(int i=0; i<m; i++) {
00209                         if(ifail[i] > 0) {
00210                                 printf("### WARNING: dsyevx_wrapper(): ifail[%d]: %d   curr_ev[%d]=%.18f\n"
00211                                         ,i, ifail[i],ifail[i],w[ifail[i]]);
00212                         }
00213                 }
00214 #endif
00215         }
00216 
00217         if ((info == 0) && (m > 0)) { // there is at least one eigenvector
00218 #ifdef TRACE_ALL
00219                 cpp_printvecDBL("eigenvalues", w, m);
00220 #endif
00221         }
00222         
00223         delete [] work;
00224         delete [] ifail;
00225         delete [] iwork;
00226 }
00227 
00228 
00229 

Generated on Wed Nov 30 03:03:59 2011 by  doxygen 1.4.7