2

有限差分法を使用して 3D シュレディンガー方程式を解くプログラムを作成しています。私のコードの 1D および 2D バージョンは問題なく実行されましたが、3D バージョンでは、行列の生成に問題があることがわかりました (QM を知っている人にとっては、これはハミルトニアン行列です。知らない人にとっては、重要ではありません) は、はるかに多くの時間を費やしています (典型的なグリッド間隔の場合は数分、最小固有値ファインダーを含む他のすべての操作の場合は数秒です!)。

私の行列生成をより効率的に記述する方法について誰か提案があるかどうか疑問に思っていました。以下に 2 つのバージョンのコードを示します。1 つは比較的理解しやすいバージョンで、もう 1 つは、スパース行列を作成するときにエントリに直接インデックスを付けるべきではなく、3 つのベクトル (行と行と列インデックス、およびそれぞれの値) を取得し、それらから疎行列を生成します。残念ながら、後者はスピードアップにはまったく役立っていません。なぜなら、私はまだばかげた三重にネストされたループを使用しており、それを回避する良い方法を思いつかないからです。

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) から 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 に答える