还是一个循环的问题

zhouds81

不好意思 再次求助！还是循环的问题，如下：

a=rand(3000,3000)
b=sparse(3*3000,3*3000)
for i=1:3000;
for j=1:3000;
      b(3*i-2,3*j-2)=a(i,j)
end
end
请问情况如何矢量化？
多谢！

想了一会，给出的解决方法如下：请各位看看是否合适
a=rand(3000,3000)；
trans=zeros(3000,3*3000);
for n=1:3000;
trans(n,3*n-2)=1
end
b=trans'*a*trans

[ 本帖最后由 zhouds81 于 2008-6-9 03:27 编辑 ]

WaitingForMe

这个问题，基本上，很难……
给你解决一下，让你体验一下什么叫矢量化。
首先的问题是：在你的代码后面增加分号，显示也是需要花时间的。
其次才是矢量化。
第一步，调试的时候，不要用这么大的矩阵来搞……
n=30;
a=rand(n,n);

tic;
b=sparse(3*n,3*n);
for i=1:n;
for j=1:n;
      b(3*i-2,3*j-2)=a(i,j);
end
end
toc

Elapsed time is 0.012402 seconds.
这个就是使用两个循环计算的时间。后面矢量化，先来了一个简单的：
tic;
c=zeros(3*n,3*n);
%b=sparse(3*n,3*n);
[am,an]=meshgrid(1:n,1:n);
[cm,cn]=meshgrid(1:3:3*n,1:3:3*n);
c(cm(:),cn(:))=a(am(:),an(:));
toc
Elapsed time is 0.116084 seconds.
判断一下结果是否正确：
>>isequal(c,b)
ans =
     1
但是矢量化以后，花费的时间反而长了很多……
所以事情都是有原因的，原因就在于，重新对稀疏矩阵赋值的时候，实际上还是需要查找的，所以结果是很慢的。
改一下就可以了，直接使用a来创建稀疏矩阵：
tic;
[dm,dn]=meshgrid(1:3:3*n,1:3:3*n);
d=sparse(dm(:),dn(:),a(:),3*n,3*n).';
toc
Elapsed time is 0.003139 seconds.
>> isequal(d,b)
ans =
     1
这下快多了！
恭喜！
最后来玩玩你这个300*300的大矩阵：
n=300;
a=rand(n,n);
tic;
[dm,dn]=meshgrid(1:3:3*n,1:3:3*n);
d=sparse(dm(:),dn(:),a(:),3*n,3*n).';
toc
Elapsed time is 0.040848 seconds.

[ 本帖最后由 bainhome 于 2008-6-9 03:58 编辑 ]

linji526

我也想了一种方法：

1. 
   
2. n=300;
   
3. a=rand(n,n);
   
4. %-----------Initial method------------
   
5. tic;
   
6. b=sparse(3*n,3*n);
   
7. for i=1:n;
   
8. for j=1:n;
   
9.       b(3*i-2,3*j-2)=a(i,j);
   
10. end
    
11. end
    
12. time0=toc
    
13. %-------------Fist method-----------
    
14. tic;
    
15. [dm,dn]=meshgrid(1:3:3*n,1:3:3*n);
    
16. d=sparse(dm(:),dn(:),a(:),3*n,3*n).';
    
17. time1=toc
    
18. %--------------Second method----------
    
19. tic;
    
20. x=1:3:3*n;
    
21. nx = length(x);
    
22. i=repmat(x,1,nx).';
    
23. j=x(ones(nx, 1),:);
    
24. j=j(:);
    
25. s=a(:);
    
26. d2=sparse(i,j,s,3*n,3*n);
    
27. tim2=toc
    
28. %-----Check equal-----
    
29. d_b=isequal(d,b)
    
30. d2_b=isequal(d2,b)
    

复制代码

三种方法比较的结果如下：
time0 =
    7.9343

time1 =
    0.0550

tim2 =
    0.0238

d_b =
     1

d2_b =
     1

zhouds81

多谢多谢。还是直接赋值比较好，尽量避免循环啊

zhouds81

多谢多谢。还是直接赋值比较好，尽量避免循环啊

rocwoods

楼主，你可真厉害，发完帖子，WaitingForMe不到一个小时就给你完整解答了，linji526也是在8小时内给你解决了，你居然再datenum('Jun-30-2010 23:41:00') - datenum('Jun-9-2008 02:50:00') =   751.8688天（两年零21.8688天）后才回应。不过比起那些别人帮助解答后再也不回的人还是要好很多。