Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

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


Math Forum » Discussions » Software » comp.soft-sys.matlab

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

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
James Tursa

Posts: 2,094
Registered: 8/5/09
Re: Lapack spptrf function
Posted: Apr 11, 2013 5:13 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

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



Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.