Fortran计算大型对称稀疏矩阵的二范数

在fortran中，有时候需要对一个大型的稀疏矩阵求取2范数。
由于大型稀疏矩阵一般使用CSR存储格式，MKL函数库中的gesvd的参数是稠密矩阵<matlab可以直接计算稀疏矩阵的2范数>
所以本文根据相关知识，给出求解CSR存储的对称方阵2范数的代码。

由矩阵分析可以知道：
特征值分解和奇异值分解的区别所有的矩阵都可以进行奇异值分解，而只有方阵才可以进行特征值分解。
当所给的矩阵是对称的方阵，A'=A，二者的结果是相同的。
也就是说对称矩阵的特征值分解是所有奇异值分解的一个特例。
但是二者还是存在一些小的差异，奇异值分解需要对奇异值从大到小的排序，而且全部是大于等于零。

代码运行环境：
win10 + vs2019 + ivf2020
代码如下：
program SparseMatrix2norm
    implicit none
    integer, parameter        :: n = 3  !// depend on your equation
    integer                   :: i, j, mm, tmp, fileid, first, num
    real(kind=8)              :: matrix(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  
    
    
    !//=====================
    character(len=1)  :: uplo = 'u'
    integer           :: fpm(128), m0, loop, m, info
    real*8, parameter :: emin = 1d-5, emax = 1d5
    real*8            :: x(n,n), epsout, e(n), res
    
    !// input data
    open( newunit = fileid, file = 'SparseMatrix.txt' )  !// The sparse matrix data needs to be given by itself
    read( fileid,* ) matrix
    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 ( *,'(<n>f7.2)' ) matrix
    write ( *,'(a)' ) '>> CSR data'
    write ( *,'(*(f7.2))' ) aa
    write ( *,'(a)' ) '>> B data'
    write ( *,'(a)' ) '>> Colnum of CSR data'
    write ( *,'(*(i4))' ) ja
    write ( *,'(a)' ) '>> The position of CSR data'
    write ( *,'(*(i4))' ) ia
  
    !// eig
    m0 = n
    call feastinit(fpm)
    call dfeast_scsrev(uplo, n, aa, ia, ja, fpm, epsout, loop, emin, emax, m0, e, x, m, res, info)
    
    if ( info == 0 ) then
            write(*,'(a)') "计算正确..."
    else
            write(*,'(a)') "计算错误..."
    end if
    
    write ( *,'(a)' ) '>> The 2-norm of Sparse Matrix is...'
    write(*,'(*(1x,f7.4,2x))') maxval(abs(e))
    
    
    deallocate( aa, ja, ia )
  
end program SparseMatrix2norm


运行结果如下：
>> Original data
   2.00   0.00  -1.00
   0.00   3.00   0.00
  -1.00   0.00   4.00
>> CSR data
   2.00  -1.00   3.00   4.00
>> B data
>> Colnum of CSR data
   1   3   2   3
>> The position of CSR data
   1   3   4   5
计算正确...
>> The 2-norm of Sparse Matrix is...
  4.4142

上述计算结果与matlab的计算结果相同。