我用matlab来实现ct算法，本质就是解线性方程组Ax=y，求出x，只是由于这个线性方程组可能是欠定的，所以使用迭代算法来求解，而代数重建算法（ART）是效果比较好的一个算法，ART的思想是：给定一个初始值，作为第一个解，将初始值代入线性方程组第一个方程，用迭代公式估算第二个解，然后将第二个估算解代入方程组第二个方程估计第三个解，以此类推，如果线性方程组有m个方程，当用第m个解代入第m个方程求出第m+1个解的估计后，则返回到第一个方程开始新一轮迭代，即用第m+1个估算解代入第一个方程求出第m+2个估算解，直到最后求出的解满足范数小于给定的误差值就退出迭代循环。ART所用的迭代公式为：
x(k+1)=x(k)+gama*a[i]'*(y[i]-a[i]*x(k))/(a[i]*a[i]')
这里，x(k)是前一次迭代的结果，x(k+1)是本次迭代的估计值，a[i]是系数矩阵A的第i行，y[i]是y列向量的第i个元素，y[i]-a[i]*x(k)表达的意义就是将前一次迭代的结果代入到线性方程组第i个方程中估算y值并与y的测量值求差，gama为一个常数，成为松弛因子。我设计的matlab代码为：

function x=ctart(A,x0,y,m,gama,tol)
% 利用ART算法解线性方程组
% CT线性方程组为 A*x=y
% 参数说明：
% A——线性方程组系数矩阵
% x0——方程组解的初始值
% y——投影数据向量
% m——投影值的个数，即向量y的维数
% gama——松弛因子
% tol——迭代容差

if (norm(y-A*x0)<=tol) % 若x0已经满足收敛条件，则退出
x=x0;
return;
end

count=1; % 计数器
ai=A(1,:); % 取A的第1行
yi=y(1,:); % 取y的第1行
beta=gama*ai'/(ai*ai'); % 构造松弛参数向量
xk=x0+beta*(yi-ai*x0) % 求第二个解

while(norm(y-A*xk)>tol) %开始迭代
x0=xk;
count=count+1; % 计数器进一位
row=mod(count,m)+1; % 计算迭代所取的方程的序号
ai=A(row,:); % 取A的第row行
yi=y(row,:); % 取y的第row行
beta=gama*ai'/(ai*ai'); % 构造松弛参数向量
xk=x0+beta*(yi-ai*x0) % 估算新的解
end

x=xk;
return;

设计完这个代码后，我代入一个四元线性方程组解算，可是解算的效率很低，计算机运算了1天半才得到结果，想请教我的代码有那些地方导致执行效率低下，可否帮我优化一下？

我是学生，现在不是很懂matlab与C++的混合编程，能指导一下么？谢谢！
我也想过用C++，所以我将上面的代码移植到了MFC上，不过很麻烦，需要创建矩阵类，代码也显得没有matlab那么简洁。我将其做成DLL然后用执行文件调用，可是运算了几个小时后就出现了错误，但又不是矩阵运算方面的错误，因为我在矩阵运算函数中都添加了错误提示，运算出错的时候程序并没有弹出我给的提示对话框，也不知道是类存泄漏问题还是什么，今晚准备话时间再运算调试一次。

[此贴子已经被作者于2007-11-13 10:06:52编辑过]

昨天又调试了一下自己编写的算法C++代码，开发工具使用的是VC，程序退出后，debug进入下面这一段
/* break into debugger at specific memory allocation */
if (lRequest == _crtBreakAlloc)
_CrtDbgBreak();
看来是和内存有关，可奇怪的是，我使用同样的矩阵类执行其它的迭代运算（比如最速下降法）并没有出现什么错误，于是我在网上搜索了一下，找到了如下文章：
http://blog.csdn.net/miracle08/archive/2006/12/24/1457060.aspx
文章大意是出现这样的错误是VC6的一个bug，如果程序（compiled with vc6 debug version run-time library）中不停地请求内存，那么该程序最终会被异常中止，只是时间问题。我的矩阵类中的成员函数就是涉及到用new来申请内存，而迭代循环中又不停地要调用矩阵类的成员函数，看来我得考虑其它方式了，现在发现VC和Matlab混编确实是个不错的选择，也在网上找了一些资料，不过，现在还比较疑惑，怎样才能做到即充分调用Matlab的数学计算，又能巧妙地让VC执行while循环呢？如果我纯粹地将算法的m文件（包括循环代码在内）全部编译成dll或者C++文件，然后使用VC调用dll的方式执行m文件内的函数，是否能够克服Matlab处理循环效率不高的弱势呢？

我用的是matlab 7.0.1，上午用共轭梯度法代码试了一下混编，好像有问题。
共轭梯度的M文件为：

function x=ctcg(A,x0,y,n,tol)
% 线性方程为y=Ax
% 参数说明：
% A——线性方程系数
% x0——初值
% y——已知数据
% n——线性方程未知数个数
% tol——迭代终止误差
% 线性方程可以化为正定二次矩阵函数f=x'*Q*x/2+b'*x+c
% 其中，Q=2*A'*A，b=-2*A'*y，c=y'*y
% 函数梯度为 gradf=Q*x+b
Q=2*A'*A;
b=-2*A'*y;
% 求初始梯度
gradf0=Q*x0+b;
k=0; % 计数器
% 若初始值已满足，则返回结果
if (norm(gradf0)<=tol)
x=x0;
return;
end
p0=-gradf0; % 初始搜索方向
lamda=-gradf0'*p0/(p0'*Q*p0); % 初始步长
xk=x0+lamda*p0; % 求第二个值
gradfk=Q*xk+b; % 求第二个值对应的梯度
k=1;
% 进入迭代运算
while (norm(gradfk)>tol)
if (k==n) % 若计数到n，则以xn为新初始值开始新一轮共轭梯度迭代
k=0; % 计数器清零
x0=xk;
gradf0=Q*x0+b;
if (norm(gradf0)<=tol)
x=x0;
return;
end
p0=-gradf0;
else
% F-R共轭梯度法
beta=norm(gradfk)*norm(gradfk)/(norm(gradf0)*norm(gradf0));
pk=-gradfk+beta*p0;
x0=xk;
gradf0=gradfk;
p0=pk;
end
lamda=-gradf0'*p0/(p0'*Q*p0); % 新的步长
xk=x0+lamda*p0; % 新的值
gradfk=Q*xk+b; % 新的梯度
k=k+1;
end
x=xk;
return;

这个代码经过线性方程组测试可以解出来解。然后mbuild和mex配置编译器，再用
mcc -B csglsharedlib:cg ctcg
来编译，这是《精通Matlab7》给的方法。编译完成后生成了dll，然后配置VC，将这个dll及其头文件和库拷贝到VC工程目录下，给一个按钮，在按钮函数中编写如下代码：

mxArray* A=NULL;
mxArray* x0=NULL;
mxArray* y=NULL;
mxArray* n=NULL;
mxArray* tol=NULL;
mxArray* x=NULL;

cgInitialize();

double da[16];
da[0] = 1.0000;
da[1] = 0.5000;
da[2] = 0.3333;
da[3] = 0.2500;
da[4] = 0.5000;
da[5] = 0.3333;
da[6] = 0.2500;
da[7] = 0.2000;
da[8] = 0.3333;
da[9] = 0.2500;
da[10] = 0.2000;
da[11] = 0.1667;
da[12] = 0.2500;
da[13] = 0.2000;
da[14] = 0.1667;
da[15] = 0.1429;

double dx0[4];
dx0[0] = 1;
dx0[1] = 1;
dx0[2] = 1;
dx0[3] = 1;

double dy[4];
dy[0] = 0.5;
dy[1] = 1;
dy[2] = 2;
dy[3] = 0;

double dn = 4;

double dtol = 0.000001;

A = mxCreateDoubleMatrix(4, 4, mxREAL);
x0 = mxCreateDoubleMatrix(4, 1, mxREAL);
y = mxCreateDoubleMatrix(4, 1, mxREAL);
n = mxCreateDoubleMatrix(1, 1, mxREAL);
tol = mxCreateDoubleMatrix(1, 1, mxREAL);
x = mxCreateDoubleMatrix(4, 1, mxREAL);

memcpy(mxGetPr(A), da, 16*sizeof(double));
memcpy(mxGetPr(x0), dx0, 4*sizeof(double));
memcpy(mxGetPr(y), dy, 4*sizeof(double));
memcpy(mxGetPr(n), &dn, sizeof(double));
memcpy(mxGetPr(tol), &dtol, sizeof(double));

mlfCtcg(1, &x, A, x0, y, n, tol);

double dx[4];
memcpy(dx, mxGetPr(x), 4*sizeof(double));

mxDestroyArray(A);
mxDestroyArray(x0);
mxDestroyArray(y);
mxDestroyArray(n);
mxDestroyArray(tol);
mxDestroyArray(x);

cgTerminate();

CString str;
str.Format("x1=%lf, x2=%lf, x3=%lf, x4=%lf", dx[0], dx[1], dx[2], dx[3]);
AfxMessageBox(str);

假设测试的线性方程组为：
x1+0.5*x2+0.3333*x3+0.25*x4=0.5
0.5*x1+0.3333*x2+0.25*x3+0.2*x4=1
0.3333*x1+0.25*x2+0.2*x3+0.1667*x4=2
0.25*x1+0.2*x2+0.1667*x3+0.1429*x4=0

最后，完成运算给出的解是：
x1=0.000000, x2=0.000000, x3=0.000000, x4=0.000000
不是我要的解。我不知道程序错在哪里……

对于第一条，在matlab 7的环境下，除了代码优化外，如何能通过混编的方式极大提升运算速度呢？
第二条，row的最大值不好说，看观测数据有多少吧，但因为方程组是欠定方程组，观测数往往没有未知数多，大概row会是即使这个量级吧，也就是说，方程组里的方程大概有即使个，但可能需要用这几十个方程去解100个左右的未知数。

PS：关于6楼的那个问题，昨天我又重新编译了一遍DLL，然后重新将生成的DLL、库和头文件添加到了VC工程中，并注意将matlab编译生成的.ctf文件也加到了工程中，最后编译VC工程，再运算，结果通过了，出来了我想要的运算结果……

[此贴子已经被作者于2007-11-15 8:53:31编辑过]

把带循环的M文件编译成DLL供VC加载是否能提高运算速度？

我在一个研究所里读书，我先试着完成，有需要帮助的时候，还得请教你。

有个问题是：我最终要发布的程序是VC编写的应用程序，现在用matlab编写算法，做成DLL，但是matlab算法里的循环又要调用mex文件，那么这么一来程序发布时还需要发布什么？将matlab的算法dll和mex的dll一起发布？