如何在Fortran中用BLAS/LAPACK由上三角生成完整对称矩阵

一、补全对称矩阵的下三角部分 #

假设你有一个N阶对称矩阵A，仅存储上三角（主对角线及以上元素有效），下三角已分配内存但未初始化。要将下三角填充为上三角的转置，得到完整对称矩阵，可通过以下BLAS操作实现：

方法1：转置+三角复制（直观易理解） #

1. 分配临时矩阵AT，维度与A一致（N×N），leading dimension设为LDA（连续存储时通常等于N）。
2. 调用BLAS的DTRAN（双精度转置）将A转置到AT：

   CALL DTRAN(N, N, A, LDA, AT, LDA)
   

   此时AT的下三角部分就是原A上三角的转置（即对称矩阵应有的下三角元素）。
3. 调用BLAS的DLACPY将AT的下三角复制到A的下三角：

   CALL DLACPY('L', N, N, AT, LDA, A, LDA)
   

   参数说明：
   • 'L'：指定复制下三角部分
   • 前两个N：矩阵的行数和列数
   • AT：源矩阵，A：目标矩阵
   • LDA：矩阵的leading dimension

方法2：直接用DSYRK（无需临时矩阵，更高效） #

如果不想额外分配临时矩阵，可直接调用DSYRK（对称秩k更新）完成填充：

CALL DSYRK('L', 'T', N, N, 1.0D0, A, LDA, 0.0D0, A, LDA)

参数说明：

• 'L'：操作目标为矩阵的下三角部分
• 'T'：使用原矩阵的转置参与计算
• 1.0D0：缩放因子alpha，0.0D0：缩放因子beta
• 该调用会将A的下三角填充为上三角的转置，最终得到完整对称矩阵。

二、修改矩阵为非对称矩阵 #

现在A已是完整对称矩阵，接下来执行两个操作：将向量x加到第一行，从第一列减去向量x。

1. 向量加到第一行 #

利用BLAS的DAXPY（双精度向量加法），由于Fortran是列主序存储，第一行元素的内存步长为LDA：

CALL DAXPY(N, 1.0D0, x, 1, A(1,1), LDA)

参数说明：

• N：向量长度
• 1.0D0：向量x的缩放因子
• x：待加的向量，1：x的元素步长（连续存储）
• A(1,1)：第一行第一个元素的起始地址，LDA：第一行元素的内存步长

2. 从第一列减去向量 #

同样用DAXPY，第一列元素是连续存储的，步长为1：

CALL DAXPY(N, -1.0D0, x, 1, A(1,1), 1)

参数说明：

• -1.0D0：向量x的缩放因子（实现减法）
• A(1,1)：第一列第一个元素的起始地址，1：第一列元素的内存步长

完整代码示例 #

PROGRAM SymToNonSym
    IMPLICIT NONE
    INTEGER, PARAMETER :: N = 5
    INTEGER :: LDA
    DOUBLE PRECISION :: A(N,N), AT(N,N), x(N)
    INTEGER :: i, j

    ! 初始化：填充上三角对称矩阵，下三角未初始化
    LDA = N
    A = 0.0D0
    DO j = 1, N
        DO i = 1, j
            A(i,j) = REAL(i + j, KIND=8)
        END DO
    END DO

    ! 初始化向量x
    x = [(REAL(i, KIND=8), i=1, N)]

    ! 步骤1：补全下三角，生成完整对称矩阵
    CALL DTRAN(N, N, A, LDA, AT, LDA)
    CALL DLACPY('L', N, N, AT, LDA, A, LDA)

    ! 步骤2：修改为非对称矩阵
    CALL DAXPY(N, 1.0D0, x, 1, A(1,1), LDA)  ! 第一行加x
    CALL DAXPY(N, -1.0D0, x, 1, A(1,1), 1)   ! 第一列减x

    ! 输出结果验证
    PRINT *, "Modified non-symmetric matrix:"
    DO i = 1, N
        PRINT *, A(i,:)
    END DO

END PROGRAM SymToNonSym

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