fortran将csr格式（全存储）的稀疏矩阵改成只存储上三角元素

chder_白南

于 2023-01-15 12:03:09 发布

阅读量492

点赞数 1

CC 4.0 BY-SA版权

分类专栏： CSR fortran 文章标签：矩阵算法线性代数

本文链接：https://blog.csdn.net/chd_lkl/article/details/128693183

fortran 同时被 2 个专栏收录

38 篇文章

订阅专栏

CSR

3 篇文章

订阅专栏

program main
        
        implicit none
        integer :: i, j, k, num, ii, jj
        !       1      2       3       4
        !-------------------------------
        !1|     10                    12
        !2|            21     33      1
        !3|            33     3            
        !4|     12     1              4
        !integer, parameter :: n = 4, nz = 10    !// n为行数, nz为非零元素个数
        !real*4  :: a(nz)   = [10.0, 12.0, 21.0, 33.0, 1.0, 33.0, 3.0, 12.0, 1.0, 4.0]
        !integer :: ia(n+1) = [1, 3, 6, 8, 11]
        !integer :: ja(nz)  = [1, 4, 2, 3, 4, 2, 3, 1, 2, 4]
        
        !// 将稀疏矩阵的数组存储上三角部分 
        !// [ 稀疏矩阵newa的大小: 7个元素 = (nz-n)/2+n ]
        !// [ 稀疏矩阵newia的大小: 5个元素 1, 3, 6, 7, 8 ]
        !       1      2       3       4
        !-------------------------------
        !1|     10                    12
        !2|            21     33      1
        !3|                   3            
        !4|                           4
        
        integer, parameter :: n = 5, nz = 13    !// n为行数, nz为非零元素个数
        real*4  :: a(nz)   = [1, 2, 4, 5, 7, 8, 9, 2, 5, 8, 11, 9, 12]
        integer :: ia(n+1) = [1, 3, 5, 8, 12, 14]
        integer :: ja(nz)  = [1, 4, 2, 4, 3, 4, 5, 1, 2, 3, 4, 3, 5]
        !       1      2       3       4        5
        !----------------------------------------
        !1|     1                      2
        !2|            4               5
        !3|                    7       8        9     
        !4|     2      5       8       11       
        !5|                    9                12
        
        !// 存储稀疏矩阵上三角部分的新数组
        real*4, allocatable  :: newa(:)
        integer              :: newnz
        integer, allocatable :: newia(:), newja(:)
        
        newnz = (nz-n)/2 + n
        allocate( newa(newnz), newja(newnz), newia(n+1) )
        
        j = 0; ii = 0; jj = 0
        newia(1) = 1
        do i = 1, n  !// 总共存储n行
                num = ia(i+1) - ia(i)   !// 原始每行的元素个数
                
                jj = 0
                do k = j+1, j+num
                        
                        if ( ja(k) >= i ) then
                                jj = jj + 1
                                ii = ii + 1
                                newja(ii) = ja(k)
                                newa(ii)  = a(k)
                        end if
                        
                end do
                newia(i+1) = newia(i) + jj
                j          = j  + num
        end do

        print*, newja
        print*, newia
        print*, newa
                        
end program main