fortran2008数组操作-4-随机矩阵

Fortran 中定义随机矩阵的方法

在 Fortran 中定义随机矩阵有多种方式，以下是几种常用的方法：

1. 使用内置随机数生成器

基本方法

program random_matrix
    implicit none
    integer, parameter :: n = 3
    real :: A(n,n)
    integer :: i, j

    ! 初始化随机数种子
    call random_seed()

    ! 生成随机矩阵
    call random_number(A)

    ! 输出矩阵
    do i = 1, n
        print '(3F10.4)', A(i,:)
    end do
end program random_matrix

指定范围的随机矩阵

program ranged_random_matrix
    implicit none
    integer, parameter :: n = 3, m = 4
    real :: A(n,m)
    real :: min_val = -1.0, max_val = 1.0

    call random_seed()
    call random_number(A)  ! 生成[0,1)区间的随机数
    A = min_val + (max_val - min_val) * A  ! 映射到[min_val, max_val)

    print *, "Random matrix in [", min_val, ",", max_val, "):"
    do i = 1, n
        print '(4F10.4)', A(i,:)
    end do
end program

2. 整数随机矩阵

program integer_random_matrix
    implicit none
    integer, parameter :: n = 5
    integer :: A(n,n)
    integer :: min_int = 0, max_int = 100

    call random_seed()
    do i = 1, n
        do j = 1, n
            call random_int(min_int, max_int, A(i,j))
        end do
    end do

    print *, "Integer random matrix:"
    do i = 1, n
        print '(5I5)', A(i,:)
    end do

contains
    subroutine random_int(lower, upper, result)
        integer, intent(in) :: lower, upper
        integer, intent(out) :: result
        real :: r

        call random_number(r)
        result = lower + floor((upper - lower + 1) * r)
    end subroutine
end program

3. 正态分布随机矩阵

使用 Box-Muller 变换

program normal_random_matrix
    implicit none
    integer, parameter :: n = 3
    real :: A(n,n)
    real :: mean = 0.0, stddev = 1.0

    call random_seed()
    do i = 1, n
        do j = 1, n
            A(i,j) = normal_random(mean, stddev)
        end do
    end do

    print *, "Normal random matrix:"
    do i = 1, n
        print '(3F10.4)', A(i,:)
    end do

contains
    function normal_random(mean, stddev) result(r)
        real, intent(in) :: mean, stddev
        real :: r
        real :: u1, u2, pi = 4.0 * atan(1.0)

        call random_number(u1)
        call random_number(u2)
        r = mean + stddev * sqrt(-2.0 * log(u1)) * cos(2.0 * pi * u2)
    end function
end program

4. 使用外部库生成随机矩阵

4.1 使用 LAPACK 的随机矩阵生成器

program lapack_random_matrix
    use, intrinsic :: iso_fortran_env, only: dp => real64
    implicit none
    integer, parameter :: n = 4
    real(dp) :: A(n,n)
    integer :: seed(4) = [42, 17, 89, 23]  ! 自定义种子

    ! 使用LAPACK的随机数生成器
    call dlarnv(3, seed, n*n, A)  ! 3表示均匀分布

    print *, "LAPACK random matrix:"
    do i = 1, n
        print '(4F10.4)', A(i,:)
    end do
end program

4.2 使用 RANDOM_NUMBER 的高性能版本

program parallel_random_matrix
    use omp_lib
    implicit none
    integer, parameter :: n = 1000
    real :: A(n,n)
    integer :: i, j

    ! 并行生成随机矩阵
    !$omp parallel do private(i,j)
    do j = 1, n
        do i = 1, n
            call random_number(A(i,j))
        end do
    end do
    !$omp end parallel do

    print *, "Generated 1000x1000 random matrix"
end program

5. 特殊类型随机矩阵

5.1 对称随机矩阵

program symmetric_random_matrix
    implicit none
    integer, parameter :: n = 4
    real :: A(n,n)
    integer :: i, j

    call random_seed()
    do i = 1, n
        do j = i, n  ! 只生成上三角部分
            call random_number(A(i,j))
            if (i /= j) A(j,i) = A(i,j)  ! 对称赋值
        end do
    end do

    print *, "Symmetric random matrix:"
    do i = 1, n
        print '(4F10.4)', A(i,:)
    end do
end program

5.2 稀疏随机矩阵

program sparse_random_matrix
    implicit none
    integer, parameter :: n = 5, sparsity = 0.3  ! 30%非零元素
    real :: A(n,n) = 0.0
    integer :: i, j
    real :: r

    call random_seed()
    do i = 1, n
        do j = 1, n
            call random_number(r)
            if (r < sparsity) call random_number(A(i,j))
        end do
    end do

    print *, "Sparse random matrix:"
    do i = 1, n
        print '(5F10.4)', A(i,:)
    end do
end program

6. 随机数种子控制

6.1 固定种子（可重复结果）

program fixed_seed_random
    implicit none
    integer, parameter :: n = 3
    real :: A(n,n)
    integer :: seed_size, clock
    integer, allocatable :: seed(:)

    ! 获取系统默认种子大小
    call random_seed(size=seed_size)
    allocate(seed(seed_size))

    ! 设置固定种子
    seed = 42  ! 可以是任意固定值
    call random_seed(put=seed)

    call random_number(A)

    print *, "Fixed seed random matrix:"
    do i = 1, n
        print '(3F10.4)', A(i,:)
    end do

    deallocate(seed)
end program

6.2 基于时间的随机种子

program time_seed_random
    implicit none
    integer, parameter :: n = 3
    real :: A(n,n)
    integer :: seed_size, clock
    integer, allocatable :: seed(:)

    ! 获取系统时间作为种子
    call system_clock(count=clock)
    call random_seed(size=seed_size)
    allocate(seed(seed_size))
    seed = clock + 37 * [(i, i=0, seed_size-1)]
    call random_seed(put=seed)

    call random_number(A)

    print *, "Time-seeded random matrix:"
    do i = 1, n
        print '(3F10.4)', A(i,:)
    end do

    deallocate(seed)
end program

7. 性能优化建议

大矩阵生成：对于大矩阵，考虑按列生成（Fortran是列优先存储）
并行生成：使用OpenMP并行化随机数生成
向量化：使用数组操作代替循环
批处理：一次性生成所有随机数比逐个生成更高效

8. 注意事项

随机数质量：内置random_number不适合加密等高质量随机需求
并行安全：某些随机数生成器在并行环境下需要特别处理
可重复性：调试时可使用固定种子
分布特性：random_number生成均匀分布，其他分布需要转换

以上方法涵盖了Fortran中定义随机矩阵的多种方式，从简单的均匀分布到更复杂的正态分布和特殊矩阵类型。根据具体需求选择最适合的方法。