代码之家 › 专栏 › 技术社区 › Jared Lo

具有多个职业的广义find()矢量化版本

vectorization gpu optimization matlab

Jared Lo · 技术社区 · 3 年前

我有下面的过程,其中N个粒子在M个格子中处理。每个盒子可以被任意数量的粒子占据。我把晶格的状态存储在一个向量中。

例如,状态向量[1 1 1 0]意味着框1中有1个粒子,框2中有1个子粒,框3中有1子粒,框4中有0个子粒。而状态向量[2 1 3 0]表示框1中的2个粒子、框2中的1个粒子、方框3中的3个粒子和框4中的0个粒子。

从这个向量中,我提取了粒子在另一个向量中的位置,以便进行进一步的计算。因此,对于状态[1 1 1 0](N<=M),位置向量为[1 2 3],这在Matlab中使用

pos=find([1 1 1 0])

对于更一般的情况向量,如[2 1 3 0](N>M),位置向量为[1 1 2 3 3]。我通过自定义的Matlab函数实现了这一点

  pos=findGeneral([2 1 3 0])

此函数的代码为

function [position]=findGeneral(state,N)

if 1*isempty(state(state>=2))==0  
    r=find(state);
    state=state(state~=0);        
    position=zeros(1,N);
    t=cell(1,length(r));    
    for i=1:length(r)
        t{i}=r(i)*ones(1,state(i));
    end
    position(1:N)=cell2mat(t);    
else
    position=find(state);    
end

end

当我的状态向量是gpuArray时,我试图运行这个函数,但显然,cell2mat()与gpuArray参数不兼容。

有没有一种方法可以对代码的7-11行进行矢量化,这样我就可以去掉for循环和cellarray了?

0 回复 | 直到 3 年前

Wolfie Radu Stefan 3 年前

您可以使用

position = [1 1 2 3 3 3];
state = sum( position.' == (1:4) );
% state = [2, 1, 3, 0]

您可以使用

state = [2, 1, 3, 0];
position = repelem( 1:numel(state), state );
% position = [1, 1, 2, 3, 3, 3]

未在gpuArray上测试(仅在普通阵列上),但从 docs ,看起来像 numel , repelem 和 sum 都由gpuArrays支持。

John Bofarull Guix 3 年前

与Wolfie的回答相反,我提供了您函数的简化版本 findGeneral 如下

1.模拟数据

ns=10
s=randi([0 3],1,ns) % s same as state
N=sum(s)

2.内容 findGeneral

y=[];
fork=1:1:ns
  s2=s(k);
  if s2==0
    y=[y 0];
  end
  while s2
    y=[y k];
    s2=s2-1;
  end
end

的回报 findGeneral 是 y 与 position 在问题中。

推荐文章

Costinean Sebastian · 如何优化使用for循环处理大型数组的MATLAB脚本的执行时间?"

1 年前

Gyro · 为什么我的%xmm3寄存器在vbroadcastsd中使用第一个参数,而不是第四个参数?

2 年前

aeiou · 对日志表进行矢量化处理,以确定最新的可用性

2 年前

IchKenneDeinenNamen · 当将向量存储为类属性时,为什么NumPy的矢量化评估速度较慢?

2 年前

roshoka · 如何在不使用for循环的情况下将函数应用于林空间的每个元素

2 年前

Xavier · 变长矢量的矢量化和

2 年前

bbbb · 矢量化或在大型数据帧上使用多重处理[重复]

2 年前

Sterling Butters · NumPy 1D近邻中循环的避免

2 年前

AnthonyML · 修改2d numpy数组中分隔符后的行中的所有元素

2 年前

Tim · np.vectorize和relativelta返回“relativeltaonlydiffs datetime/date”

2 年前