From: Markus <mb_matlabREMOVE@freenetTHIS.de>
Path: news.mathworks.com!newsfeed-00.mathworks.com!webcrossing
Newsgroups: comp.soft-sys.matlab
Subject: Re: DSYEVX functionality in MATLAB
Message-ID: <ef5e907.0@webcrossing.raydaftYaTP>
Date: Tue, 24 Jul 2007 07:22:39 -0400
References: <f83o4n$fth$1@login.nmt.edu>
Lines: 65
NNTP-Posting-Host: 129.69.32.222
MIME-Version: 1.0
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
Xref: news.mathworks.com comp.soft-sys.matlab:420696



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;
}