lapack

Easily call any LAPACK or BLAS routine from inside Matlab.

Description

This package contains the command "lapack", which provides a simple interface to call any LAPACK or BLAS routine from inside Matlab, as well as the command "lapackhelp", which brings up detailed information about any LAPACK or BLAS routine.

A summary of the features include:
· Uses Matlab's internal LAPACK and BLAS libraries.
· Automatic argument type conversion (ie. converts arguments from Matlab types to FORTRAN types).
· Easy access to detailed descriptions about the behavior and arguments for every LAPACK and BLAS function
· Allows use of LAPACK and BLAS without any knowledge of MEX, C, or FORTRAN.

Acknowledgements

This file inspired Trilin and Tridiegs.

• MATLAB release: MATLAB 7.9 (R2009b)

[dm]

 

[Dhrubo]

Also, support for PLASMA routines would be nice.

[Dhrubo]

Hi, I have a request. Could you please also add support for SCALAPACK routines?

Regards,
Dhrubo

[Felipe]

Qianyong Chen, I think the segfault was because lwork was too small; doing as follows doesn't seem to crash matlab. The fortranError as well as fortranError_SVD remain both numerically zero.

A = rand(118,106);

% 1st: results from direct matlab
[UM, SM, VM] = svd(A);

matlabError = max(max(abs(A-UM*SM*VM')))

% 2nd: results when using 'lapack' by Tim Toolan
[m,n] = size(A);

C = lapack('DGESVD', 'A','A',m, ...
n,A,m, zeros(n,1), zeros(m),m,zeros(n),n,zeros(1,1),-1,0)
lwork = C{end-2}(1);
C = lapack('DGESVD', 'A','A',m, ...
n,A,m, zeros(n,1), zeros(m),m,zeros(n),n,zeros(lwork,1),lwork,0);
[SF, UF, VFt] = C{[7,8,10]};
VF = VFt';

fortranError_U = max(max(abs(UM - UF)))
fortranError_V = max(max(abs(VM - VF)))
fortranError_SVD = max(abs(SF - diag(SM) ))

fortranError = max(max(abs( A - UF * SM * VF')))

[Dhrubo]

Dude, your code rocks! You're awesome! :) I have been able to use the dsyev routine on a very large matrix (size 200k^2). It actually worked just as well, if not better, than Matlab 2011b's eig.

[Qianyong Chen]

I really like this idea. But couldn't make it work on my machine: 64bit Win7 with SDK 7.1 with Matlab 2011b.
The first run is always fine. But when I run the same code a few more times, it either stops due to segmentation error,
or gives unreliable results. Also, the difference for U and V (between Matlab results and using lapack.c provided here)
is large when the matrix is as small as 118 by 106. The error for the matrix itself (A - U * S * V) is very small though.

What I did is: 1) compile with mex lapack.c
2) run the following code:

% test a couple subroutines
clear all;
%--------------------
%1st: dgesvd.f
%--------------------
A = rand(118,106);

% 1st: results from direct matlab
[UM, SM, VM] = svd(A);

matlabError = max(max(abs(A-UM*SM*VM')))

% 2nd: results when using 'lapack' by Tim Toolan
[m,n] = size(A);

C = lapack('DGESVD', 'A','A',m, ...
n,A,m, zeros(n,1), zeros(m),m,zeros(n),n,zeros(5*m,1),5*m,0);
[SF, UF, VF] = C{[7,8,10]};

fortranError_U = max(max(abs(UM - UF)))
fortranError_V = max(max(abs(VM - VF')))
fortranError_SVD = max(abs(SF - diag(SM) ))

fortranError = max(max(abs( A - UF * SM * VF)))

[Felipe G. Nievinski]

 

[Felipe G. Nievinski]

Tim, I'm using a slightly modified version of your lapack package in my trilin package (#31125). I'd be glad if you'd consider incorporating the modifications in my submission into your original submission. Thanks.

[Felipe G. Nievinski]

I was suuuuper excited about this package. I am considering ditching all my LAPACK-calling C mex code in favor of it. The compilation is a bliss, no hassle in finding out where either or both LAPACK and BLAS are installed, depending on MATLAB version, computer, and compiler!

Before that, though, I was wondering if could you please comment on the following issues.

- Compiling without the -largeArrayDims flag in a 64-bit MATLAB seems to give unreliable results. Any idea why? It should be optional.

- The instructions on top of lapack.c say:
To build on Windows, run: mex lapack.c
Since lapack.m already beautifully takes care of that using largeArrayDims, you might want to change lapack.c to instruct users to simply run lapack(.m).

- Using lapack.mexw64 compiled on MATLAB 2009 64-bit gives unreliable results when used on MATLAB 2007 64-bit; it seems related to casting integer types, e.g., ipiv in a call to dgetrf has some very large values:
17179869189
21474836483
0
0

- Empty input crashes MATLAB, e.g., lapack('dgetrf', 0, 0, [], 0, 0, 1)

- I see that all arguments are copied to new memory. For me that's a concern, because I need to handle very large matrices. I understand doing that for arguments that will be overwritten by LAPACK, but that's clearly unnecessary for, e.g., the matrix factor in solving a linear system. I was wondering if you'd agree with the inclusion of a new argument, something like this:

else if (c == 'U'){
/* USER DEFINED **********************************************
/ User guarantees type and also whether or not LAPACK will overwrite this argument.
/ lowercase: input-only: need not copy.
/ UPPERCASE: input/output: need to copy.
/ TODO: complex-valued variables need at least to have their real/imag parts reorganized. */
if (mxIsComplex(args[i])){
sprintf(msg, "Can't convert argument %d to FORTRAN type "
"DOUBLE PRECISION, use double.", i+1);
mexErrMsgTxt(msg);
}
if (isupper(argType[i])){
value = mxDuplicateArray(args[i]);
a[i] = mxGetPr(value);
if (plhs[0] && mxIsCell(plhs[0]))
mxSetCell(plhs[0], indArgOut[i], value);
else
plhs[indArgOut[i]] = value;
}
else{
a[i] = mxGetPr(args[i]);
}
}

- e.g., lapackhelp('dtrcon') leads to an unreachable website; the most stable URL seems to be:
http://www.netlib.org/lapack/double/dtrcon.f
(where the double/single subdir would have to be chosen automatically)

Thanks for your time.

[Anastasios]

Thanks for sharing your work with us, Tim.

Best,
Tasos

[worker b]

I'm trying to edit my previous review. I am now able to get this to run! To compile on 2006b under Gentoo linux:
mex -D_GNU_SOURCE -ldl -lmwlapack lapack.c
I have heard that -lmwlapack is not needed for >2006 versions of matlab, tho.
I look forward to using this! (the initial timings look good, too.)

[worker b]

Unable to get this to run on my machine, but this is probably due to a general matlab/lapack issue. However, I had to use the -D_GNU_SOURCE flag in compiling (RTLD_DEFAULT is a GNU extension evidently. whatever.)