单进程初始化全局矩阵后分发至进程网格调用ScaLAPACK的内存问题

问题根源 #

1. 内存耗尽原因：原代码中所有进程均分配了全局矩阵mat(nst,nst)，当nst>20000时，单进程内存占用约3.2GB（双精度浮点），多进程叠加后直接超出可用RAM限制。
2. 初始化错误原因：仅根进程初始化mat时，其他进程的mat未定义，后续通过indxl2g读取mat(i,j)会获取未初始化的垃圾值，导致本地子矩阵Aloc赋值错误，最终PDSYEV计算结果异常。

最优解决方案：本地直接计算子矩阵元素（完全消除全局矩阵） #

核心思路是让每个进程仅计算自己负责的本地子矩阵Aloc的元素，无需分配全局矩阵，从根源避免内存浪费，同时保证所有进程的Aloc正确初始化。

修正后代码

! 假设已完成BLACS进程网格初始化（nprow, npcol, myrow, mycol, ictxt 已定义）
block = 64  ! 根据CPU缓存调整，推荐64/128
lld = numroc(nst, block, myrow, 0, nprow)  ! 计算本地矩阵首维度

! 初始化ScaLAPACK描述符
call descinit(descAloc, nst, nst, block, block, 0, 0, ictxt, lld, info)
call descinit(descz, nst, nst, block, block, 0, 0, ictxt, lld, info)

! 计算工作区大小（保留原逻辑）
lworkk = max(1 + 6*nst + 2*np*nq, 3*nst + max(block*(np+1), 3*block)) + 2*nst

! 仅分配本地所需数组，**删除全局矩阵mat和Zglobal**
allocate(Aloc(lld,lld), workk(lworkk), wrk(np), Z(lld,lld))

! 每个进程直接计算本地子矩阵元素，无需全局矩阵
do local_i = 1, np
    ! 本地索引转全局索引
    i = indxl2g(local_i, block, myrow, 0, nprow)
    do local_j = 1, nq
        j = indxl2g(local_j, block, mycol, 0, npcol)
        ! 直接计算全局(i,j)位置的元素值
        if (i == j) then
            Aloc(local_i, local_j) = real(rpsi(i)**2/3.0 + rpsi(i)**3/2.0)
        else
            Aloc(local_i, local_j) = real(rpsi(i)**3/2.0)
        end if
    end do
end do

! 调用ScaLAPACK特征值求解
call pdsyev("V", "U", nst, Aloc, 1, 1, descAloc, W, Z, 1, 1, descz, workk, lworkk, info)
write(0,*) "info pdsyev", info 

! 根进程输出结果
if (myrow == 0 .and. mycol == 0) then 
    write(0,*) "eigenvalues in parallel"
    write(0,"(6f8.2)") W(1), W(2), W(3), W(nst-2), W(nst-1), W(nst) 
end if

! 释放内存
deallocate(Aloc, workk, wrk, Z)

方案优势

• 内存占用最小：每个进程仅持有本地子矩阵，内存占用降至lld*lld*8字节（lld通常远小于nst）。
• 效率更高：避免了进程间数据分发的开销，直接本地计算。
• 无初始化错误：所有进程的Aloc均通过计算生成，无垃圾值。

备选方案：根进程初始化全局矩阵后分发 #

若无法直接本地计算元素（如元素依赖全局矩阵的其他位置），可采用根进程初始化全局矩阵，再分发到各进程的方式：

! 仅根进程分配并初始化全局矩阵
if (myrow == 0 .and. mycol == 0) then
    allocate(mat(nst,nst))
    mat = 0.d0
    do i=1,nst
        mat(i,i) = mat(i,i) + rpsi(i)**2/3.0
        do j=1,nst
            mat(i,j) = mat(i,j) + rpsi(i)**3/2.0
        end do
    end do
    ! 初始化全局矩阵的ScaLAPACK描述符
    call descinit(descMat, nst, nst, nst, nst, 0, 0, ictxt, nst, info)
end if

! 分配本地数组（不含全局矩阵）
allocate(Aloc(lld,lld), workk(lworkk), wrk(np), Z(lld,lld))

! 根进程将全局矩阵分发到各进程的Aloc
call pdgemr2d(nst, nst, mat, 1, 1, descMat, Aloc, 1, 1, descAloc, ictxt, info)

! 根进程释放全局矩阵
if (myrow == 0 .and. mycol == 0) then
    deallocate(mat)
end if

! 后续PDSYEV调用与输出同最优方案

关键注意事项 #

• 确保BLACS进程网格正确初始化（需调用blacs_pinfo、blacs_gridinit等函数）。
• 块大小block需匹配CPU缓存，通常选64或128，平衡计算与通信开销。
• 本地矩阵首维度lld需通过numroc函数计算，避免数组越界。
• 所有ScaLAPACK描述符必须通过descinit正确初始化，参数错误会导致计算异常。

内容的提问来源于stack exchange，提问作者FortranZT