53 lines
1.4 KiB
C
53 lines
1.4 KiB
C
|
/************************************************************************
|
||
|
|
* U = mexchol(X) X = U'*U
|
||
|
|
* Computes upper-triangular Cholesky factor U of full symmetric positive
|
||
|
|
* definite matrix X.
|
||
|
|
*
|
||
|
|
* Uses LAPACK routine DPOTRF
|
||
|
|
************************************************************************/
|
||
|
|
|
||
|
|
#include "mex.h"
|
||
|
|
#include <math.h>
|
||
|
|
|
||
|
|
/********************************************************************
|
||
|
|
PROCEDURE mexFunction - Entry for Matlab
|
||
|
|
*********************************************************************/
|
||
|
|
void mexFunction(const int nlhs, mxArray *plhs[],
|
||
|
|
const int nrhs, const mxArray *prhs[])
|
||
|
|
{
|
||
|
|
int n;
|
||
|
|
double *x, *info;
|
||
|
|
char *chn;
|
||
|
|
int infotmp=0, j,k,jn;
|
||
|
|
|
||
|
|
|
||
|
|
if(nrhs < 1)
|
||
|
|
mexErrMsgTxt("mexchol: requires 1 input arguments.");
|
||
|
|
if(nlhs > 1)
|
||
|
|
mexErrMsgTxt("mexchol: requires 1 output argument.");
|
||
|
|
|
||
|
|
x = mxGetPr(prhs[0]);
|
||
|
|
if (mxIsSparse(prhs[0]))
|
||
|
|
mexErrMsgTxt("Sparse X not supported by mexchol.");
|
||
|
|
n = mxGetM(prhs[0]);
|
||
|
|
if ( n != mxGetN(prhs[0]) )
|
||
|
|
mexErrMsgTxt("X should be square.");
|
||
|
|
plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
|
||
|
|
info = mxGetPr(plhs[0]);
|
||
|
|
|
||
|
|
chn = "U";
|
||
|
|
dpotrf(chn, &n, x, &n, &infotmp);
|
||
|
|
info[0] = (double)infotmp;
|
||
|
|
/***** make lower triangular part of x zero *****/
|
||
|
|
for (j=0; j<n; j++) {
|
||
|
|
jn = j*n;
|
||
|
|
for (k=j+1; k<n; k++) { x[k+jn] = 0; }
|
||
|
|
}
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
/************************************************************************/
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|