"James Tursa" wrote in message <kk5ur5$gml$1@newscl01ah.mathworks.com>...
> "Francesco Perrone" <francesco86perrone@yahoo.it> wrote in message <kk5ner$r67$1@newscl01ah.mathworks.com>...
> > Dear Mr. Tursa,
> >
> > I have a certain n*n matrix. Then calling
> >
> > my_matrix = nn_matrix(itril(size(nn_matrix)));
> >
> > I am able to wrap it into a lower triangular factor.
> >
> > Afterwards, I manipulate my_matrix by means of a coherence function; the outcome, for only one step of the above mentioned FOR loop, is
> >
> > Coh_u = exp(a.*my_matrix.*sqrt((f(ii)/Uhub).^2 + (0.12/Lc).^2)).*(df.*psd(ii,1));
> >
> > Therefore, Coh_u is alread packed as spptrf input.
> >
> > Then I call the lapack wrapper as follows:
> >
> > [C1u,C2u,HH_u,C3u] = lapack('spptrf','L',size(nn_matrix,1),Coh_u,0);
> >
> > HH_u is the expected result.
> >
> > I would more than appreciate if you could write a mex file of the spptrf function for me. Unfortunately I started only now learning some C++ and it would require a long time for me to catch up with mex file.
> >
> > Look really forward to hearing from you.
>
> Here is the mex function. It assumes you are passing it a lower triangular matrix in packed form (either single or double) and returns the same. You will have to link in the BLAS and LAPACK libraries. Here is how you would do it on 32bit Windows:
>
> liblapack = [matlabroot '\extern\lib\win32\microsoft\libmwlapack.lib'];
> libblas = [matlabroot '\extern\lib\win32\microsoft\libmwblas.lib'];
> mex('spptrf.c',liblapack,libblas)
>
> or on earlier versions of MATLAB that do not have a separate BLAS library:
>
> liblapack = [matlabroot '\extern\lib\win32\microsoft\libmwlapack.lib'];
> mex('spptrf.c',liblapack)
>
> If you are using the lcc compiler instead of a Microsoft compiler then use lcc instead of microsoft in the directory strings above.
>
> On Unix systems you will have to link in the libraries differently. Consult your doc for these instructions.
>
> James Tursa
>
> 
>
> /* File: spptrf.c
> * Syntax: B = spptrf(A)
> * [B,INFO] = spptrf(A)
> * Purpose: Computes the LAPACK function SPPTRF (or DPPTRF) of the input, which
> * is assumed to be in lower triangular packed form on input and output
> * Programmer: James Tursa
> * Date: 10Apr2013
> */
>
> #include "mex.h"
>
> #if defined(__OS2__)  defined(__WINDOWS__)  defined(WIN32)  defined(_WIN32)  defined(WIN64)  defined(_WIN64)  defined(_MSC_VER)
> #define SPPTRF spptrf
> #define DPPTRF dpptrf
> #else
> #define SPPTRF spptrf_
> #define DPPTRF dpptrf_
> #endif
>
> #ifndef MWSIZE_MAX
> #define mwSignedIndex int
> #endif
>
> void SPPTRF(char *, mwSignedIndex *, float *, mwSignedIndex *);
> void DPPTRF(char *, mwSignedIndex *, double *, mwSignedIndex *);
>
> void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
> {
> char UPLO = 'L';
> mwSignedIndex INFO, M, N = 1;
>
> if( nlhs > 2 ) {
> mexErrMsgTxt("Too many outputs.");
> }
> if( nrhs != 1 
> mxGetNumberOfDimensions(prhs[0]) != 2 
> (mxGetM(prhs[0]) != 1 && mxGetN(prhs[0]) != 1) 
> (!mxIsSingle(prhs[0]) && !mxIsDouble(prhs[0])) 
> mxIsComplex(prhs[0])  mxIsEmpty(prhs[0])  mxIsSparse(prhs[0]) ) {
> mexErrMsgTxt("Need exactly one full real single or double vector input.");
> }
> plhs[0] = mxDuplicateArray(prhs[0]);
> M = mxGetNumberOfElements(prhs[0]);
> while( (N*(N+1))/2 < M ) { /* lazy way to do this */
> N++;
> }
> if( (N*(N+1))/2 != M ) {
> mexErrMsgTxt("Input vector is not the proper length.");
> }
> if( mxIsSingle(prhs[0]) ) {
> SPPTRF( &UPLO, &N, mxGetData(plhs[0]), &INFO );
> } else {
> DPPTRF( &UPLO, &N, mxGetData(plhs[0]), &INFO );
> }
> if( nlhs == 2 ) {
> plhs[1] = mxCreateDoubleScalar(INFO);
> }
> }
Dear Mr. Tursa,
I really thank for writing the mex wrapper for spptrf. I will test as soon as I will have access to MATLAB.
Kindest regards,
Francesco Perrone
