2

我正在编写一个程序来使用有限差分法求解 3D Schroedinger 方程。我的代码的 1D 和 2D 版本运行得很好，但是在 3D 版本中，我发现矩阵的生成（对于那些知道 QM 的人，这是哈密顿矩阵；对于那些不知道的人，这并不重要）到目前为止花费的时间最多（典型的网格间距需要几分钟，而所有其他操作需要几秒钟，包括最小特征值查找器！）。

我想知道是否有人对如何更有效地编写矩阵生成有任何建议。我在下面包含了我的代码的两个版本：一个应该相对容易理解，然后是遵循 MATLAB 文档建议的第二个版本，我不应该在制作稀疏矩阵时直接索引条目，而是应该制作三个向量（行和列索引及其各自的值）并从中生成稀疏矩阵。不幸的是，后者根本没有帮助加快速度，因为我仍在使用愚蠢的三重嵌套循环，而且我想不出避免它的好方法。

delta = 0.1e-9;
Lx = 2e-9;
x = 0:delta:Lx;
Nx = length(x);
Ly = 2e-9;
y = 0:delta:Ly;
Ny = length(y);
Lz = 2e-9;
z = 0:delta:Lz;
Nz = length(z);

map = inline('((idx_x-1) * Ny*Nz) + ((idx_y-1) * Nz) + idx_z','idx_x','idx_y','idx_z','Ny','Nz'); % define an inline helper function for mapping (x,y,z) indices to a linear index

Tsparse = sparse([],[],[],Nx*Ny*Nz, Nx*Ny*Nz, 7*(Nx-2)*(Ny-2)*(Nz-2)); % kinetic part of Hamiltonian matrix: (d^2/dx^2 + d^2/dy^2 + d^2/dz^2); NOTE: we'll have 7*(Nx-2)*(Ny-2)*(Nz-2) non-zero entries in this matrix, so we get the sparse() function to preallocate enough memory for this

for idx_x = 2:Nx-1
    for idx_y = 2:Ny-1
        for idx_z = 2:Nz-1
            Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x ,idx_y , idx_z , Ny, Nz) ) = -6/delta^2;
            Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x+1,idx_y , idx_z , Ny, Nz) ) = 1/delta^2;
            Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x-1,idx_y , idx_z , Ny, Nz) ) = 1/delta^2;
            Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x ,idx_y+1, idx_z , Ny, Nz) ) = 1/delta^2;
            Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x ,idx_y-1, idx_z , Ny, Nz) ) = 1/delta^2;
            Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x ,idx_y , idx_z+1, Ny, Nz) ) = 1/delta^2;
            Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x ,idx_y , idx_z-1, Ny, Nz) ) = 1/delta^2;
        end
   end
end

此代码创建了一个矩阵，该矩阵仅沿 7 个对角线具有非零条目（并且并非每个对角线中的所有条目都非零）。

这是我尝试以更接近 MATLAB 文档建议的方式创建 T 矩阵的代码版本：

delta = 0.1e-9;
Lx = 2e-9;
x = 0:delta:Lx;
Nx = length(x);
Ly = 2e-9;
y = 0:delta:Ly;
Ny = length(y);
Lz = 2e-9;
z = 0:delta:Lz;
Nz = length(z);

map = inline('((idx_x-1) * Ny*Nz) + ((idx_y-1) * Nz) + idx_z','idx_x','idx_y','idx_z','Ny','Nz'); % define an inline helper function for mapping (x,y,z) indices to a linear index

Iidx = zeros(7*(Nx-2)*(Ny-2)*(Nz-2),1); % matrix row indices
Jidx = zeros(7*(Nx-2)*(Ny-2)*(Nz-2),1); % matrix col indices
vals = zeros(7*(Nx-2)*(Ny-2)*(Nz-2),1); % matrix non-zero values, corresponding to (row,col)
cnt = 1;
for idx_x = 2:Nx-1
    for idx_y = 2:Ny-1
        for idx_z = 2:Nz-1
            % Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x ,idx_y , idx_z , Ny, Nz) ) = -6/delta^2;
            Iidx(cnt) = map(idx_x,idx_y,idx_z,Ny,Nz);
            Jidx(cnt) = map(idx_x,idx_y,idx_z,Ny,Nz);
            vals(cnt) = -6/delta^2;
            cnt = cnt + 1;
            % Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x+1,idx_y , idx_z , Ny, Nz) ) = 1/delta^2;
            Iidx(cnt) = map(idx_x,idx_y,idx_z,Ny,Nz);
            Jidx(cnt) = map(idx_x+1,idx_y,idx_z,Ny,Nz);
            vals(cnt) = 1/delta^2;
            cnt = cnt + 1;
            % Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x-1,idx_y , idx_z , Ny, Nz) ) = 1/delta^2;
            Iidx(cnt) = map(idx_x,idx_y,idx_z,Ny,Nz);
            Jidx(cnt) = map(idx_x-1,idx_y,idx_z,Ny,Nz);
            vals(cnt) = 1/delta^2;
            cnt = cnt + 1;
            % Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x ,idx_y+1, idx_z , Ny, Nz) ) = 1/delta^2;
            Iidx(cnt) = map(idx_x,idx_y,idx_z,Ny,Nz);
            Jidx(cnt) = map(idx_x,idx_y+1,idx_z,Ny,Nz);
            vals(cnt) = 1/delta^2;
            cnt = cnt + 1;
            % Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x ,idx_y-1, idx_z , Ny, Nz) ) = 1/delta^2;
            Iidx(cnt) = map(idx_x,idx_y,idx_z,Ny,Nz);
            Jidx(cnt) = map(idx_x,idx_y-1,idx_z,Ny,Nz);
            vals(cnt) = 1/delta^2;
            cnt = cnt + 1;
            % Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x ,idx_y , idx_z+1, Ny, Nz) ) = 1/delta^2;
            Iidx(cnt) = map(idx_x,idx_y,idx_z,Ny,Nz);
            Jidx(cnt) = map(idx_x,idx_y,idx_z+1,Ny,Nz);
            vals(cnt) = 1/delta^2;
            cnt = cnt + 1;
            % Tsparse( map(idx_x,idx_y,idx_z,Ny,Nz) , map(idx_x ,idx_y , idx_z-1, Ny, Nz) ) = 1/delta^2;
            Iidx(cnt) = map(idx_x,idx_y,idx_z,Ny,Nz);
            Jidx(cnt) = map(idx_x,idx_y,idx_z-1,Ny,Nz);
            vals(cnt) = 1/delta^2;
            cnt = cnt + 1;

        end
    end
end
Tsparse = sparse(Iidx, Jidx, vals, Nx*Ny*Nz, Nx*Ny*Nz);

在此先感谢您的任何建议！

--dx.dy.dz

（旁注：“map”函数用于从 3D 坐标系 (x,y,z) 到 1D 值。假设我的特征值问题是 H psi = E psi，其中 H 是哈密顿矩阵，而 psi是一个向量，而 E 是一个标量。矩阵 H = T + V（V 未在代码示例中显示，只有 T 是）写在一个基础上，其中 3D psi 函数被离散化并从 3D 折叠到 1D。例如，假设我每个维度只有 2 个网格点，所以 x=1:1:2, y=1:​​1:2, z=1:1:2。然后我的哈密顿量写在基础 {psi( 1,1,1), psi(1,1,2), psi(1,2,1), psi(1,2,2), psi(2,1,1), psi(2,1,2) ), psi(2,2,1), psi(2,2,2)}，即它是一个 8×8 矩阵。eigs() 求解器输出的特征向量 psi 将是一个 8 分量向量，如果需要，我可以将其重新整形为 2x2x2 矩阵。）

4

1 回答 1

1

我想我可以给一些建议：

  • 而不是您自己的地图，您可能会考虑该sub2ind功能

  • 您反复调用map(idx_x,idx_y,idx_z,Ny,Nz)-确保您可以将其存储起来以供重复使用。

  • 邻居的相对位置也将保持不变 - 无需重新计算

作为一个小例子，我会这样做：

siz = [4,4,4];

pos = sub2ind(siz,1,1,1)

tmp = [
    sub2ind(siz,2,1,1)-pos
    sub2ind(siz,1,2,1)-pos
    sub2ind(siz,1,1,2)-pos
    ];

neighbors = [tmp;-tmp];
%%
big_dim = prod(siz);
mat = sparse(big_dim,big_dim);
%%
for i=2:siz(1)-1
    for j=2:siz(2)-1
        for k=2:siz(3)-1
            c_pos = sub2ind(siz,i,j,k);
            mat(c_pos,c_pos)=-6;
            c_neighbors=c_pos+neighbors;
            mat(c_pos,c_neighbors)=1;
        end
    end
end
于 2013-01-18T12:34:07.300 回答