Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: Lapack spptrf function
Replies: 10   Last Post: Apr 12, 2013 11:08 AM

 Messages: [ Previous | Next ]
 James Tursa Posts: 2,326 Registered: 8/5/09
Re: Lapack spptrf function
Posted: Apr 11, 2013 5:13 AM

"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 32-bit 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: 10-Apr-2013
*/

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

Date Subject Author
4/9/13 Francesco Perrone
4/10/13 Steven Lord
4/10/13 Francesco Perrone
4/10/13 Steven Lord
4/10/13 James Tursa
4/10/13 James Tursa
4/11/13 Francesco Perrone
4/11/13 James Tursa
4/11/13 Francesco Perrone
4/12/13 Francesco Perrone
4/12/13 James Tursa