在有限差分求解器中高效地生成稀疏矩陣

我正在編寫一個程序來使用有限差分方法求解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

（附註：「地圖」功能用於從3D去座標系（X，Y，Z），以一維值假設我的特徵值問題H psi = E psi，其中H是哈密爾頓矩陣，而psi是矢量，E是標量。矩陣H = T + V（代碼樣本中未示出V，只有T是）被寫入例如，假設每個維度只有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的矩陣，一個特徵向量psi eigs（）求解器的輸出將是一個8分量矢量，然後我就可以解析了。高原肺水腫回，如果我想要一個2x2x2的矩陣）

來源

2013-01-18 dxdydz

我想我可以給幾個指針：

而不是你自己的地圖，你可能會考慮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-18 12:34:07 bdecaf

在有限差分求解器中高效地生成稀疏矩陣

回答

相關問題