欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

fortran恢复CSR格式的数组

程序员文章站 2023-12-25 22:24:21
...
program RestoreCSRdata
    implicit none
    integer, parameter        :: n = 5
    integer                   :: i, j, k, mm, tmp, fileid, first, num
    real(kind=8)              :: matrix(n,n), x(n), a(n,n)
    real(kind=8), allocatable :: aa(:)
    integer, allocatable      :: ia(:), ja(:)
  
    type                      :: node
        real(kind=8)          :: s
        integer               :: k1, k2
        type (node), pointer  :: next
    end type node 
    type (node), pointer      :: head, tail, p, q  
  
    !// input data
    open ( newunit = fileid, file = 'SparseMatrix.txt' )  !// The sparse matrix data needs to be given by itself
    do i = 1, n
            read ( fileid,* ) matrix(i,:)
    end do
    close ( fileid )

    !// =========================store data by CSR format=======================
    first = 1
    if ( .not.associated(p) ) then
        allocate( p )
        q    => p
        nullify( p%next )
        q%k2 = first
        tmp  = q%k2
    end if
  
    num = 0  !// calculate the number of no-zero
    do i =  1, n
        mm = 0  !// calculate the position of no-zero in every row
        do j = i, n
            if ( matrix(i,j) /= 0.0 ) then
                if ( .not.associated(head) ) then
                    allocate( head )
                    tail    => head
                    nullify( tail%next )
                    tail%s  = matrix(i,j)
                    tail%k1 = j
                    num     = num + 1
                    mm      = mm + 1
                else
                    allocate( tail%next )
                    tail    => tail%next
                    tail%s  = matrix(i,j)
                    tail%k1 = j
                    num     = num + 1  
                    mm      = mm + 1  
                    nullify( tail%next )
                end if
            end if
        end do
        if ( associated(p) ) then
            allocate( q%next )
            q    => q%next
            q%k2 = tmp + mm
            tmp  = q%k2
            nullify( q%next )
        end if
    end do
  
    allocate( aa(num), ja(num), ia(n+1) )  !// store the no-zero element
    !// output a and ja
    tail => head
    j = 1
    do 
        if ( .not.associated(tail) ) exit
        aa(j) = tail%s
        ja(j) = tail%k1
        j     = j + 1
        tail  => tail%next
    end do
  
    !// output ia
    q => p
    j = 1
    do 
        if ( .not.associated(q) ) exit
        ia(j) = q%k2
        j     = j + 1
        q     => q%next
    end do
  
    nullify( head, tail, p, q )
    !// =========================store data by CSR format=======================
    write ( *,'(a)' ) '>> Original data'
    write ( *,'(5f7.2)' ) matrix
    write ( *,'(a)' ) '>> CSR data'
    write ( *,'(*(f7.2))' ) aa
    write ( *,'(a)' ) '>> Colnum of CSR data'
    write ( *,'(*(i4))' ) ja
    write ( *,'(a)' ) '>> The position of CSR data'
    write ( *,'(*(i4))' ) ia
  
  
    !// Restore the original array
    a = 0.d0
    j = 0
    do i = 1, n
                    num = ia(i+1) - ia(i)
                    do k = j+1, j+num
                        a(i,ja(k)) = aa(k)
                        a(ja(k),i) = aa(k)
                    end do
                    j = j + num
    end do
    deallocate( aa, ja, ia )
    
    write ( *,'(a)' ) '>> Restore the original array'
    do i = 1, n
            write(*,'(*(f7.2))') a(i,:)
    end do
  
end program RestoreCSRdata

运行结果如下:
>> Original data
   1.00   0.00   0.00   2.00   0.00
   0.00   4.00   0.00   5.00   0.00
   0.00   0.00   7.00   8.00   9.00
   2.00   5.00   8.00  11.00   0.00
   0.00   0.00   9.00   0.00  12.00
>> CSR data
   1.00   2.00   4.00   5.00   7.00   8.00   9.00  11.00  12.00
>> Colnum of CSR data
   1   4   2   4   3   4   5   4   5
>> The position of CSR data
   1   3   5   8   9  10
>> Restore the original array
   1.00   0.00   0.00   2.00   0.00
   0.00   4.00   0.00   5.00   0.00
   0.00   0.00   7.00   8.00   9.00
   2.00   5.00   8.00  11.00   0.00
   0.00   0.00   9.00   0.00  12.00
请按任意键继续. . .

 

相关标签: fortran

上一篇:

下一篇: