Thread Subject: DSYEVX functionality in MATLAB

The eig() command provides easy access to the functionality of the
LAPACK DSYEV routine. However, there are times when it is useful to
use options available only through the DSYEVX expert driver. For
example, if I'm interested only in those eigenvalues of a matrix
in a particular interval, DSYEVX can find them.
 
I had hoped that eig() would have options to get at this functionality,
but it appears that it doesn't. Is there some relatively easy way to
do this in MATLAB? I realize that this could be done with a mex file,
but that seems like a lot of work...

--
Brian Borchers borchers@nmt.edu
Department of Mathematics http://www.nmt.edu/~borchers/
New Mexico Tech Phone: 505-835-5813
Socorro, NM 87801 FAX: 505-835-5366

Using a mex-function is a bit of work, but this is definitely a
possible and fast way. Below I have pasted a mex-function that calls
the Lapack function DPOTRF which computes a Cholesky factorization
(i.e. imitating the functionality of Matlab-function chol). This not
exactly what you want, but with some copy and paste and the
documentation of your function on <http://www.netlib.org/lapack/>,
ou should quickly get your wanted function.

Just put the code below into a function (e.g. called cholmex.c),
start "mex -setup", select the gcc compiler under Linux or Lcc under
Windows and run "mex chomex.c" to compile your mex-function - ready.

If you need more help, just let me know.

Markus

#include <mex.h>
#include <matrix.h>

#ifndef WIN32
extern void dpotrf_(char *uplo, int *n, double *a, int *lda, int
*info);
#define dpotrf dpotrf_
#endif

void chol(double *matrix, int N);

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray
*prhs[] )
{
int N;

/* Input arguments */
N = mxGetN(prhs[0]);

if(N != mxGetM(prhs[0]))
mexErrMsgTxt("Matrix must be square.");

/* Output arguments */
plhs[0] = mxDuplicateArray(prhs[0]);

/* Cholesky factorization using LAPACK function */
chol(mxGetPr(plhs[0]), N);
}

void chol(double *matrix, int N)
{
  char uplo = 'U';
  int info, row, col;
  
/* Cholesky factorization using LAPACK function */
dpotrf(&uplo, &N, matrix, &N, &info);

/* Check result of function call */
if(info != 0)
{
mexPrintf("Return code of function DPOTRF: %d\n", info);
mexErrMsgTxt("Error using LAPACK function.");
}

/* delete lower triangle elements */
for(row=1; row<N; row++)
for(col=0; col<row; col++)
matrix[row + N*col] = 0;
}

Using a mex-function is a bit of work, but this is definitely a
possible and fast way. Below I have pasted a mex-function that calls
the Lapack function DPOTRF which computes a Cholesky factorization
(i.e. imitating the functionality of Matlab-function chol). This not
exactly what you want, but with some copy and paste and the
documentation of your function on <http://www.netlib.org/lapack/>,
ou should quickly get your wanted function.

Just put the code below into a function (e.g. called cholmex.c),
start "mex -setup", select the gcc compiler under Linux or Lcc under
Windows and run "mex chomex.c" to compile your mex-function - ready.

If you need more help, just let me know.

Markus

#include <mex.h>
#include <matrix.h>

#ifndef WIN32
extern void dpotrf_(char *uplo, int *n, double *a, int *lda, int
*info);
#define dpotrf dpotrf_
#endif

void chol(double *matrix, int N);

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray
*prhs[] )
{
int N;

/* Input arguments */
N = mxGetN(prhs[0]);

if(N != mxGetM(prhs[0]))
mexErrMsgTxt("Matrix must be square.");

/* Output arguments */
plhs[0] = mxDuplicateArray(prhs[0]);

/* Cholesky factorization using LAPACK function */
chol(mxGetPr(plhs[0]), N);
}

void chol(double *matrix, int N)
{
  char uplo = 'U';
  int info, row, col;
  
/* Cholesky factorization using LAPACK function */
dpotrf(&uplo, &N, matrix, &N, &info);

/* Check result of function call */
if(info != 0)
{
mexPrintf("Return code of function DPOTRF: %d\n", info);
mexErrMsgTxt("Error using LAPACK function.");
}

/* delete lower triangle elements */
for(row=1; row<N; row++)
for(col=0; col<row; col++)
matrix[row + N*col] = 0;
}

/*
If you test this function, note that DPOTRF
expects a real, symmetric, positive definite
matrix as input. For example use the following
line: x = eye(5)+rand(5); cholmex(x+x');
/*

Sorry for accidentally posting two times. Refer to my SECOND reply,
it contains some additional text.

Markus

Did anybody try to translate the mex-file above under Windows? There you have to give the path to the lapack library additionally:

mex('cholmex.c', [matlabroot, '/extern/lib/win32/lcc/libmwlapack.lib'])

If your file is called "cholmex.c" and you are using the LCC compiler.

Markus