http://www.mathworks.com/matlabcentral/newsreader/view_thread/328184
MATLAB Central Newsreader  Lapack spptrf function
Feed for thread: Lapack spptrf function
enus
©19942015 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Tue, 09 Apr 2013 14:07:20 +0000
Lapack spptrf function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328184#901934
Francesco Perrone
Hi everybody,<br>
<br>
I have got a question concerning a specific LAPACK function: spptrf.<br>
<br>
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.<br>
<br>
Therefore is the given matrix has N*N dimension, spptrf will receive as input only the N(N+1)/2 independent element stored as unidimensional matrix.<br>
<br>
Now, I've seen that C:\Program Files\MATLAB\R2012a\extern\include\lapack.h includes spptrf. But can I directly call spptrf?<br>
<br>
I do need this function in order to drastically reduce the computational effort of my current code.<br>
<br>
I thank you in advance for any support you will show.

Wed, 10 Apr 2013 15:24:43 +0000
Re: Lapack spptrf function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328184#902032
Steven_Lord
<br>
<br>
"Francesco Perrone" <francesco86perrone@yahoo.it> wrote in message <br>
news:kk17ao$o6s$1@newscl01ah.mathworks.com...<br>
> Hi everybody,<br>
><br>
> I have got a question concerning a specific LAPACK function: spptrf.<br>
><br>
> spptrf allows for a Cholesky factorization of a given positive definite <br>
> symmetric matrix by only input of the upper (respectively lower) <br>
> triangular factor, stored as an array.<br>
><br>
> Therefore is the given matrix has N*N dimension, spptrf will receive as <br>
> input only the N(N+1)/2 independent element stored as unidimensional <br>
> matrix.<br>
><br>
> Now, I've seen that C:\Program Files\MATLAB\R2012a\extern\include\lapack.h <br>
> includes spptrf. But can I directly call spptrf?<br>
<br>
Why not simply call CHOL?<br>
<br>
You could probably write a MEXfile to do this, though there are a few <br>
gotchas that you should review in the MEX documentation. But again, if <br>
you're looking for the Cholesky decomposition, there's a function in MATLAB <br>
for that called CHOL.<br>
<br>
 <br>
Steve Lord<br>
slord@mathworks.com<br>
To contact Technical Support use the Contact Us link on <br>
<a href="http://www.mathworks.com">http://www.mathworks.com</a>

Wed, 10 Apr 2013 17:42:13 +0000
Re: Lapack spptrf function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328184#902052
Francesco Perrone
Dear Mr. Lord,<br>
<br>
yes I am aware of the function CHOL and I am already using it.<br>
<br>
But CHOL gets tremendously slow when factorizing, e.g., a 1945*1945 matrix.<br>
<br>
Within my code I call CHOL, within a 2048 iterations FOR loop, for three different matrices with the above mentioned size.<br>
<br>
The elapsed time is about an hour.<br>
<br>
For the same case, after downloading the LAPACK wrapper (<a href="http://www.mathworks.com/matlabcentral/fileexchange/16777lapack">http://www.mathworks.com/matlabcentral/fileexchange/16777lapack</a>), the elapsed time was around 17 min less than before: quite a bargain I'd urge.<br>
<br>
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.

Wed, 10 Apr 2013 19:47:40 +0000
Re: Lapack spptrf function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328184#902070
Steven_Lord
<br>
<br>
"Francesco Perrone" <francesco86perrone@yahoo.it> wrote in message <br>
news:kk489l$fnr$1@newscl01ah.mathworks.com...<br>
> Dear Mr. Lord,<br>
><br>
> yes I am aware of the function CHOL and I am already using it.<br>
><br>
> But CHOL gets tremendously slow when factorizing, e.g., a 1945*1945 <br>
> matrix.<br>
<br>
Please define "tremendously slow." When I use a sample SPD matrix on my <br>
machine, it takes about a tenth of a second. YTMV.<br>
<br>
<br>
A = gallery('minij', 1945);<br>
tic<br>
[R, p] = chol(A);<br>
toc<br>
<br>
<br>
> Within my code I call CHOL, within a 2048 iterations FOR loop, for three <br>
> different matrices with the above mentioned size.<br>
><br>
> The elapsed time is about an hour.<br>
<br>
So 3*2048 = 6144 calls to CHOL, averaging about half a second per call, <br>
assuming that's the only thing in the loop? Depending on the hardware <br>
configuration of your machine, that may not be too bad.<br>
<br>
> For the same case, after downloading the LAPACK wrapper <br>
> (<a href="http://www.mathworks.com/matlabcentral/fileexchange/16777lapack">http://www.mathworks.com/matlabcentral/fileexchange/16777lapack</a>), the <br>
> elapsed time was around 17 min less than before: quite a bargain I'd urge.<br>
><br>
> Btw, after a bunch of tests, MATLAB started crashing at any time I tried <br>
> calling the function: blatantly, I don't even have a clue why that <br>
> happened.<br>
<br>
By "the function" do you mean CHOL or the LAPACK wrapper function from the <br>
File Exchange? If the former, please send the crash log, a MATfile <br>
containing the matrix, and the exact command you used to Technical Support <br>
so we can investigate the cause of the crash. If you mean the LAPACK wrapper <br>
function, that suggests that it's not operating correctly. I'd be cautious <br>
about using it until you and/or the developer can determine the cause of the <br>
crash and correct it.<br>
<br>
 <br>
Steve Lord<br>
slord@mathworks.com<br>
To contact Technical Support use the Contact Us link on <br>
<a href="http://www.mathworks.com">http://www.mathworks.com</a><br>
YTMV = 'Your Timing May Vary"

Wed, 10 Apr 2013 23:19:11 +0000
Re: Lapack spptrf function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328184#902086
James Tursa
"Francesco Perrone" <francesco86perrone@yahoo.it> wrote in message <kk489l$fnr$1@newscl01ah.mathworks.com>...<br>
> Dear Mr. Lord,<br>
> <br>
> yes I am aware of the function CHOL and I am already using it.<br>
> <br>
> But CHOL gets tremendously slow when factorizing, e.g., a 1945*1945 matrix.<br>
> <br>
> Within my code I call CHOL, within a 2048 iterations FOR loop, for three different matrices with the above mentioned size.<br>
> <br>
> The elapsed time is about an hour.<br>
> <br>
> For the same case, after downloading the LAPACK wrapper (<a href="http://www.mathworks.com/matlabcentral/fileexchange/16777lapack">http://www.mathworks.com/matlabcentral/fileexchange/16777lapack</a>), the elapsed time was around 17 min less than before: quite a bargain I'd urge.<br>
> <br>
> 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.<br>
<br>
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 copyincopyout which could be improved upon with a custom mex routine.<br>
<br>
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.<br>
<br>
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.<br>
<br>
James Tursa

Wed, 10 Apr 2013 23:59:11 +0000
Re: Lapack spptrf function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328184#902089
James Tursa
"James Tursa" wrote in message <kk4s1f$g38$1@newscl01ah.mathworks.com>...<br>
> <br>
> 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.<br>
<br>
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.<br>
<br>
James Tursa

Thu, 11 Apr 2013 07:07:07 +0000
Re: Lapack spptrf function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328184#902109
Francesco Perrone
Dear Mr. Tursa,<br>
<br>
I have a certain n*n matrix. Then calling <br>
<br>
my_matrix = nn_matrix(itril(size(nn_matrix)));<br>
<br>
I am able to wrap it into a lower triangular factor.<br>
<br>
Afterwards, I manipulate my_matrix by means of a coherence function; the outcome, for only one step of the above mentioned FOR loop, is<br>
<br>
Coh_u = exp(a.*my_matrix.*sqrt((f(ii)/Uhub).^2 + (0.12/Lc).^2)).*(df.*psd(ii,1));<br>
<br>
Therefore, Coh_u is alread packed as spptrf input.<br>
<br>
Then I call the lapack wrapper as follows:<br>
<br>
[C1u,C2u,HH_u,C3u] = lapack('spptrf','L',size(nn_matrix,1),Coh_u,0);<br>
<br>
HH_u is the expected result.<br>
<br>
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.<br>
<br>
Look really forward to hearing from you.

Thu, 11 Apr 2013 09:13:09 +0000
Re: Lapack spptrf function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328184#902120
James Tursa
"Francesco Perrone" <francesco86perrone@yahoo.it> wrote in message <kk5ner$r67$1@newscl01ah.mathworks.com>...<br>
> Dear Mr. Tursa,<br>
> <br>
> I have a certain n*n matrix. Then calling <br>
> <br>
> my_matrix = nn_matrix(itril(size(nn_matrix)));<br>
> <br>
> I am able to wrap it into a lower triangular factor.<br>
> <br>
> Afterwards, I manipulate my_matrix by means of a coherence function; the outcome, for only one step of the above mentioned FOR loop, is<br>
> <br>
> Coh_u = exp(a.*my_matrix.*sqrt((f(ii)/Uhub).^2 + (0.12/Lc).^2)).*(df.*psd(ii,1));<br>
> <br>
> Therefore, Coh_u is alread packed as spptrf input.<br>
> <br>
> Then I call the lapack wrapper as follows:<br>
> <br>
> [C1u,C2u,HH_u,C3u] = lapack('spptrf','L',size(nn_matrix,1),Coh_u,0);<br>
> <br>
> HH_u is the expected result.<br>
> <br>
> 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.<br>
> <br>
> Look really forward to hearing from you.<br>
<br>
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:<br>
<br>
liblapack = [matlabroot '\extern\lib\win32\microsoft\libmwlapack.lib'];<br>
libblas = [matlabroot '\extern\lib\win32\microsoft\libmwblas.lib'];<br>
mex('spptrf.c',liblapack,libblas)<br>
<br>
or on earlier versions of MATLAB that do not have a separate BLAS library:<br>
<br>
liblapack = [matlabroot '\extern\lib\win32\microsoft\libmwlapack.lib'];<br>
mex('spptrf.c',liblapack)<br>
<br>
If you are using the lcc compiler instead of a Microsoft compiler then use lcc instead of microsoft in the directory strings above.<br>
<br>
On Unix systems you will have to link in the libraries differently. Consult your doc for these instructions.<br>
<br>
James Tursa<br>
<br>
<br>
<br>
/* File: spptrf.c<br>
* Syntax: B = spptrf(A)<br>
* [B,INFO] = spptrf(A)<br>
* Purpose: Computes the LAPACK function SPPTRF (or DPPTRF) of the input, which<br>
* is assumed to be in lower triangular packed form on input and output<br>
* Programmer: James Tursa<br>
* Date: 10Apr2013<br>
*/<br>
<br>
#include "mex.h"<br>
<br>
#if defined(__OS2__)  defined(__WINDOWS__)  defined(WIN32)  defined(_WIN32)  defined(WIN64)  defined(_WIN64)  defined(_MSC_VER)<br>
#define SPPTRF spptrf<br>
#define DPPTRF dpptrf<br>
#else<br>
#define SPPTRF spptrf_<br>
#define DPPTRF dpptrf_<br>
#endif<br>
<br>
#ifndef MWSIZE_MAX<br>
#define mwSignedIndex int<br>
#endif<br>
<br>
void SPPTRF(char *, mwSignedIndex *, float *, mwSignedIndex *);<br>
void DPPTRF(char *, mwSignedIndex *, double *, mwSignedIndex *);<br>
<br>
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])<br>
{<br>
char UPLO = 'L';<br>
mwSignedIndex INFO, M, N = 1;<br>
<br>
if( nlhs > 2 ) {<br>
mexErrMsgTxt("Too many outputs.");<br>
}<br>
if( nrhs != 1 <br>
mxGetNumberOfDimensions(prhs[0]) != 2 <br>
(mxGetM(prhs[0]) != 1 && mxGetN(prhs[0]) != 1) <br>
(!mxIsSingle(prhs[0]) && !mxIsDouble(prhs[0])) <br>
mxIsComplex(prhs[0])  mxIsEmpty(prhs[0])  mxIsSparse(prhs[0]) ) {<br>
mexErrMsgTxt("Need exactly one full real single or double vector input.");<br>
}<br>
plhs[0] = mxDuplicateArray(prhs[0]);<br>
M = mxGetNumberOfElements(prhs[0]);<br>
while( (N*(N+1))/2 < M ) { /* lazy way to do this */<br>
N++;<br>
}<br>
if( (N*(N+1))/2 != M ) {<br>
mexErrMsgTxt("Input vector is not the proper length.");<br>
}<br>
if( mxIsSingle(prhs[0]) ) {<br>
SPPTRF( &UPLO, &N, mxGetData(plhs[0]), &INFO );<br>
} else {<br>
DPPTRF( &UPLO, &N, mxGetData(plhs[0]), &INFO );<br>
}<br>
if( nlhs == 2 ) {<br>
plhs[1] = mxCreateDoubleScalar(INFO);<br>
}<br>
}

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

Fri, 12 Apr 2013 12:54:06 +0000
Re: Lapack spptrf function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328184#902241
Francesco Perrone
Dear Mr. Tursa,<br>
<br>
I just made a little change to your code:<br>
<br>
between lin 44 and 46 your called a while loop:<br>
<br>
while( (N*(N+1))/2 < M ) { /* lazy way to do this */<br>
N++;<br>
}<br>
and I replaced it with a do...while one:<br>
<br>
do{<br>
N++;<br>
} while ( N*(N+1)/2 < M );<br>
<br>
Do you think it will increase anyhow the computation performances?

Fri, 12 Apr 2013 15:08:08 +0000
Re: Lapack spptrf function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328184#902255
James Tursa
"Francesco Perrone" <francesco86perrone@yahoo.it> wrote in message <kk905e$ose$1@newscl01ah.mathworks.com>...<br>
> Dear Mr. Tursa,<br>
> <br>
> I just made a little change to your code:<br>
> <br>
> between lin 44 and 46 your called a while loop:<br>
> <br>
> while( (N*(N+1))/2 < M ) { /* lazy way to do this */<br>
> N++;<br>
> }<br>
> and I replaced it with a do...while one:<br>
> <br>
> do{<br>
> N++;<br>
> } while ( N*(N+1)/2 < M );<br>
> <br>
> Do you think it will increase anyhow the computation performances?<br>
<br>
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.<br>
<br>
James Tursa