Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Lapack spptrf function

Subject: Lapack spptrf function

From: Francesco Perrone

Date: 9 Apr, 2013 14:07:20

Message: 1 of 11

Hi everybody,

I have got a question concerning a specific LAPACK function: spptrf.

spptrf allows for a Cholesky factorization of a given positive definite symmetric matrix by only input of the upper (respectively lower) triangular factor, stored as an array.

Therefore is the given matrix has N*N dimension, spptrf will receive as input only the N(N+1)/2 independent element stored as uni-dimensional matrix.

Now, I've seen that C:\Program Files\MATLAB\R2012a\extern\include\lapack.h includes spptrf. But can I directly call spptrf?

I do need this function in order to drastically reduce the computational effort of my current code.

I thank you in advance for any support you will show.

Subject: Lapack spptrf function

From: Steven_Lord

Date: 10 Apr, 2013 15:24:43

Message: 2 of 11



"Francesco Perrone" <francesco86perrone@yahoo.it> wrote in message
news:kk17ao$o6s$1@newscl01ah.mathworks.com...
> Hi everybody,
>
> I have got a question concerning a specific LAPACK function: spptrf.
>
> spptrf allows for a Cholesky factorization of a given positive definite
> symmetric matrix by only input of the upper (respectively lower)
> triangular factor, stored as an array.
>
> Therefore is the given matrix has N*N dimension, spptrf will receive as
> input only the N(N+1)/2 independent element stored as uni-dimensional
> matrix.
>
> Now, I've seen that C:\Program Files\MATLAB\R2012a\extern\include\lapack.h
> includes spptrf. But can I directly call spptrf?

Why not simply call CHOL?

You could probably write a MEX-file to do this, though there are a few
gotchas that you should review in the MEX documentation. But again, if
you're looking for the Cholesky decomposition, there's a function in MATLAB
for that called CHOL.

--
Steve Lord
slord@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Subject: Lapack spptrf function

From: Francesco Perrone

Date: 10 Apr, 2013 17:42:13

Message: 3 of 11

Dear Mr. Lord,

yes I am aware of the function CHOL and I am already using it.

But CHOL gets tremendously slow when factorizing, e.g., a 1945*1945 matrix.

Within my code I call CHOL, within a 2048 iterations FOR loop, for three different matrices with the above mentioned size.

The elapsed time is about an hour.

For the same case, after downloading the LAPACK wrapper (http://www.mathworks.com/matlabcentral/fileexchange/16777-lapack), the elapsed time was around 17 min less than before: quite a bargain I'd urge.

Btw, after a bunch of tests, MATLAB started crashing at any time I tried calling the function: blatantly, I don't even have a clue why that happened.

Subject: Lapack spptrf function

From: Steven_Lord

Date: 10 Apr, 2013 19:47:40

Message: 4 of 11



"Francesco Perrone" <francesco86perrone@yahoo.it> wrote in message
news:kk489l$fnr$1@newscl01ah.mathworks.com...
> Dear Mr. Lord,
>
> yes I am aware of the function CHOL and I am already using it.
>
> But CHOL gets tremendously slow when factorizing, e.g., a 1945*1945
> matrix.

Please define "tremendously slow." When I use a sample SPD matrix on my
machine, it takes about a tenth of a second. YTMV.


A = gallery('minij', 1945);
tic
[R, p] = chol(A);
toc


> Within my code I call CHOL, within a 2048 iterations FOR loop, for three
> different matrices with the above mentioned size.
>
> The elapsed time is about an hour.

So 3*2048 = 6144 calls to CHOL, averaging about half a second per call,
assuming that's the only thing in the loop? Depending on the hardware
configuration of your machine, that may not be too bad.

> For the same case, after downloading the LAPACK wrapper
> (http://www.mathworks.com/matlabcentral/fileexchange/16777-lapack), the
> elapsed time was around 17 min less than before: quite a bargain I'd urge.
>
> Btw, after a bunch of tests, MATLAB started crashing at any time I tried
> calling the function: blatantly, I don't even have a clue why that
> happened.

By "the function" do you mean CHOL or the LAPACK wrapper function from the
File Exchange? If the former, please send the crash log, a MAT-file
containing the matrix, and the exact command you used to Technical Support
so we can investigate the cause of the crash. If you mean the LAPACK wrapper
function, that suggests that it's not operating correctly. I'd be cautious
about using it until you and/or the developer can determine the cause of the
crash and correct it.

--
Steve Lord
slord@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com
YTMV = 'Your Timing May Vary"

Subject: Lapack spptrf function

From: James Tursa

Date: 10 Apr, 2013 23:19:11

Message: 5 of 11

"Francesco Perrone" <francesco86perrone@yahoo.it> wrote in message <kk489l$fnr$1@newscl01ah.mathworks.com>...
> Dear Mr. Lord,
>
> yes I am aware of the function CHOL and I am already using it.
>
> But CHOL gets tremendously slow when factorizing, e.g., a 1945*1945 matrix.
>
> Within my code I call CHOL, within a 2048 iterations FOR loop, for three different matrices with the above mentioned size.
>
> The elapsed time is about an hour.
>
> For the same case, after downloading the LAPACK wrapper (http://www.mathworks.com/matlabcentral/fileexchange/16777-lapack), the elapsed time was around 17 min less than before: quite a bargain I'd urge.
>
> Btw, after a bunch of tests, MATLAB started crashing at any time I tried calling the function: blatantly, I don't even have a clue why that happened.

1) Which exact function(s) did you call with the lapack FEX package? And what syntax did you use to call it? I would need to double check this, but my recollection is that this package does a copy-in-copy-out which could be improved upon with a custom mex routine.

2) If you called spptrf, did you do the packing and unpacking on the MATLAB side before you called it? And did you pass it a single class variable? If you intend to pass a double class variable you need to use dpptrf instead.

3) If you need help setting up a custom mex routine to do this I can probably write it for you pretty quickly. Just ask if you need this.

James Tursa

Subject: Lapack spptrf function

From: James Tursa

Date: 10 Apr, 2013 23:59:11

Message: 6 of 11

"James Tursa" wrote in message <kk4s1f$g38$1@newscl01ah.mathworks.com>...
>
> 3) If you need help setting up a custom mex routine to do this I can probably write it for you pretty quickly. Just ask if you need this.

I just did a quick test with a custom mex routine. Calling dpptrf is taking 50% longer than calling chol. My guess is this is reflecting the necessary packing and unpacking of the data in order to call it (all that data copying takes time). So unless you can keep all the data in a packed state through all your iterations I doubt a custom mex routine using dpptrf (or spptrf as the case may be) will help your timing.

James Tursa

Subject: Lapack spptrf function

From: Francesco Perrone

Date: 11 Apr, 2013 07:07:07

Message: 7 of 11

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.

Subject: Lapack spptrf function

From: James Tursa

Date: 11 Apr, 2013 09:13:09

Message: 8 of 11

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

Subject: Lapack spptrf function

From: Francesco Perrone

Date: 11 Apr, 2013 12:42:07

Message: 9 of 11

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

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

Subject: Lapack spptrf function

From: Francesco Perrone

Date: 12 Apr, 2013 12:54:06

Message: 10 of 11

Dear Mr. Tursa,

I just made a little change to your code:

between lin 44 and 46 your called a while loop:

while( (N*(N+1))/2 < M ) { /* lazy way to do this */
 N++;
 }
and I replaced it with a do...while one:

do{
    N++;
    } while ( N*(N+1)/2 < M );

Do you think it will increase anyhow the computation performances?

Subject: Lapack spptrf function

From: James Tursa

Date: 12 Apr, 2013 15:08:08

Message: 11 of 11

"Francesco Perrone" <francesco86perrone@yahoo.it> wrote in message <kk905e$ose$1@newscl01ah.mathworks.com>...
> Dear Mr. Tursa,
>
> I just made a little change to your code:
>
> between lin 44 and 46 your called a while loop:
>
> while( (N*(N+1))/2 < M ) { /* lazy way to do this */
> N++;
> }
> and I replaced it with a do...while one:
>
> do{
> N++;
> } while ( N*(N+1)/2 < M );
>
> Do you think it will increase anyhow the computation performances?

That loop is likely in the noise as far as performance goes. The only thing you are doing with your modification is skipping the N=1 case, which you will not even be able to measure in the performance. The only reason I used a loop in the first place is because I didn't feel like using a sqrt to extract N. You could always pass in N in the argument list as the 2nd input to avoid this loop altogether if you wanted.

James Tursa

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us