在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的计算结果相同。