最近在利用matlab计算的时候遇到了一些精度上的问题。考虑到matlab本身是调用LAPACK库做矩阵运算的，所以我也想直接调用LAPACK里一些专门的库函数来处理一些特殊矩阵的运算，解决这个问题。
需要做的一些准备：
1. 阅读matlab帮助里关于C-interface的调用部分
2. 下载LAPACK库的source code
3. www.netlib.org/lapack上有lapack的user guide
matlab是解释型语言，它的运算速度和经过编译的c，fortran程序相比要慢的多，但相对的好处是，它的直观性决定了它是一种很好的建模工具，你可以用它搭起一个基本的框架，然后把需要speed up的portion用C,Fortran或者J***A来rewrite,并通过interface和原先的matlab程序联系起来。这样，你可以更加容易的在实现基本功能的基础上对程序作进一步的改进。毕竟，一个好的程序的第一步，是它能用，而不是一个代码写的多漂亮的绣花枕头，用matlab写程序就是基于这样的基本想法吧。
好了，言归正传。在matlab里提供了一个基本的interface接口，这里我只讨论c程序接口。在matlab中用mex来编译C接口程序，这样，生成的*.mexglx文件可以作为matlab的内建函数直接调用。在$matlab/extern/examples下面提供了许多有用的调用例程，一开始，可以通过阅读它们获得对c-interface文件格式的感性认识，比如，它提供了自己的指针类型，如何实现和fortran里数据存储格式的转化等等，需要特别小心。下面，我仅以两个例程来作为说明：
1：对一个复厄米矩阵求逆。
代码如下：
＃i nclude "mex.h"
＃i nclude "matrix.h"
＃i nclude
double* mat2fort(const mxArray *X, int ldz, int ndz);
mxArray* fort2mat(double *Z, int ldz, int m,int n);
void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
/*这里nlhs是输出的变量个数，plhs[]是输出的变量，nrhs:输入变量个数，prhs[]输入变量*/
{
char UPLO;
long N, *IPIV, LWORK, INFO;
double *A, *WORK;
UPLO = 'L';
N = mxGetN(prhs[0]);
IPIV = mxCalloc(N, sizeof(long));
LWORK = N*N;
A = mat2fort(prhs[0], N, N);
WORK = mxCalloc(2*N*N, sizeof(double));/*分配空间应该考虑复数的空间分配*/
/* Call complex LAPACK function */
/*关于函数部分应该参考LAPACK source code里具体的说明和user guide*/
zhetrf_ ( &UPLO, &N, A, &N, IPIV, WORK, &LWORK, &INFO);
zhetri_ ( &UPLO, &N, A, &N, IPIV, WORK, &INFO);
/* Convert output to MATLAB format */
plhs[0] = fort2mat(A, N, N, N);
/* Free intermediate Fortran format arrays */
mxFree(A); mxFree(IPIV); mxFree(WORK);
}
/*
* Convert MATLAB complex matrix to Fortran complex storage.
* Z = mat2fort(X,ldz,ndz) converts MATLAB's mxArray X to Fortran's
* complex*16 Z(ldz,ndz). The parameters ldz and ndz determine the
* storage allocated for Z, while mxGetM(X) and mxGetN(X) determine
% the amount of data copied.
*/
double* mat2fort(
const mxArray *X,
int ldz,
int ndz
)
{
int i,j,m,n,incz,cmplxflag;
double *Z,*xr,*xi,*zp;
Z = (double *) mxCalloc(2*ldz*ndz, sizeof(double));
xr = mxGetPr(X);
xi = mxGetPi(X);
m = mxGetM(X);
n = mxGetN(X);
zp = Z;
incz = 2*(ldz-m);
cmplxflag = (xi != NULL);
for (j = 0; j < n; j++) {
if (cmplxflag) {
for (i = 0; i < m; i++) {
*zp++ = *xr++;
*zp++ = *xi++;
}
} else {
for (i = 0; i < m; i++) {
*zp++ = *xr++;
zp++;
}
}
zp += incz;
}
return(Z);
}
/*
* Convert Fortran complex storage to MATLAB real and imaginary parts.
* X = fort2mat(Z,ldz,m,n) copies Z to X, producing a complex mxArray
* with mxGetM(X) = m and mxGetN(X) = n.
*/
mxArray* fort2mat(
double *Z,
int ldz,
int m,
int n
)
{
int i,j,incz;
double *xr,*xi,*zp;
mxArray *X;
X = mxCreateDoubleMatrix(m,n,mxCOMPLEX);
xr = mxGetPr(X);
xi = mxGetPi(X);
zp = Z;
incz = 2*(ldz-m);
for (j = 0; j < n; j++) {
for (i = 0; i < m; i++) {
*xr++ = *zp++;
*xi++ = *zp++;
}
zp += incz;
}
return(X);
}
好，让我们把这个c文件叫inv1.c。用mex编译，我们就可以在matlab里象一个函数那样调用它了。
这是matlab里得到的运行结果:
>> A=rand(3)+rand(3)*i;A=A+A'
A =
1.0392 0.9742 - 0.4878i 0.4850 - 0.0119i
0.9742 + 0.4878i 1.1776 0.6960 - 0.4768i
0.4850 + 0.0119i 0.6960 + 0.4768i 0.4331
>> inv(A)
ans =
0.3669 - 0.0000i 0.1432 - 0.7901i 0.2288 + 1.4373i
0.1432 + 0.7901i -0.3907 - 0.0000i 0.4457 - 1.3109i
0.2288 - 1.4373i 0.4457 + 1.3109i -0.0670 - 0.0000i
>> C=tril(inv1(A));C+C'-diag(diag(C))%这里我们用的是下三角LU分解来求逆，所以要再作一点变换
ans =
0.3669 0.1432 - 0.7901i 0.2288 + 1.4373i
0.1432 + 0.7901i -0.3907 0.4457 - 1.3109i
0.2288 - 1.4373i 0.4457 + 1.3109i -0.0670
我们用C-interface直接调用LAPACK计算得到了和matlab同样的结果。
我们对一个500x500的矩阵做10次计算，这是测试的时间比较：
matlab cputime: 9seconds
self lapack cputime: 6seconds
可以看到，效率明显提高了。
甚至，如果要做进一步提高效率，我们甚至可以用ifc做优化，重新编译LAPACK，调用优化后的库函数，当然，这是后续的工作。。。
同样的，我们还可以用LAPACK库来处理本征值问题，方法类似。但使用LAPACK库是一件挺辛苦的事，尤其是要做复杂的事，需要处理的选项就多了。。。。所以，一般只是对特殊的问题我们需要这样做，比如，如果我们处理的是一般的矩阵，matlab本身的效率并不算太低，在我的测试下，和调用 LAPACK里general矩阵的函数的效果差不多。所以如果不是特别有必要，象我这样的懒人，也就还是拿matlab来说事了，呵呵。
最后，感谢亦斌提供的有益信息，让我少走了不少弯路，呼呼。
最新评论