在fortran中，有時候需要對乙個大型的稀疏矩陣求取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 ( *,'(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的計算結果相同。
對稱矩陣 稀疏矩陣的壓縮儲存
對稱矩陣 稀疏矩陣的壓縮儲存 1 對稱矩陣的壓縮儲存 對稱矩陣顧名思義就是符合行和列的個數相同，並且矩陣中儲存的資料上三角和下三角中對應位置上的元素值是相等的。為了能夠減少儲存的空間，我們可以只儲存上三角矩陣 或者下三角矩陣中的元素，這樣就能夠極大地節省空間的浪費。下面是對稱矩陣的示列 假設對稱矩陣...
資料結構 稀疏矩陣，對稱矩陣
稀疏矩陣 在矩陣中，若數值為0的元素數目遠遠多於非0元素的數目時，則稱該矩陣為稀疏矩陣 include include include using namespace std templatestruct triple templateclass sparsematrix sparsematrix ...
矩陣 對稱矩陣及其壓縮儲存 稀疏矩陣
什麼是對稱矩陣 symmetricmatrix 對稱對稱 看 設乙個n n的方陣a，a中任意元素aij，當且僅當aij aji 0 i n 1 0 j n 1 則矩陣a是對稱矩陣。以矩陣的對角線為分隔，分為上三角和下三角。壓縮存就是矩陣儲存時只需要儲存上三角 下三角的資料，所以最多儲存n n 1 2...