fortran MKL関数ライブラリを用いて大型疎行列の特徴値と特徴ベクトルを計算する
5149 ワード
program scsrev
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), tempx(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)
tempx = x
if ( info == 0 ) then
write(*,'(a)') " ..."
else
write(*,'(a)') " ..."
end if
write ( *,'(a)' ) '>> eigenvalue is...'
write(*,'(*(1x,f7.4,2x))') e
write ( *,'(a)' ) '>> eigenvector is...'
do i = 1, n
write(*,'(*(1x,f7.4,2x))') x(i,:)
end do
do i = 1, n !// V*inv(D)
x(:,i) = x(:,i) / e(i)
end do
write ( *,'(a)' ) '>> '
x = matmul(x, transpose(tempx))
do i = 1, n
write(*,'(*(1x,f7.4,2x))') x(i,:) !// x
end do
!//
write ( *,'(a)' ) '>> '
x = matmul(matrix, x)
do i = 1, n
write(*,'(*(1x,f7.4,2x))') x(i,:) !//
end do
deallocate( aa, ja, ia )
end program scsrev
:
>> 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
...
>> eigenvalue is...
1.5858 3.0000 4.4142
>> eigenvector is...
-0.9239 0.0000 -0.3827
-0.0000 1.0000 0.0000
-0.3827 0.0000 0.9239
>>
0.5714 0.0000 0.1429
0.0000 0.3333 0.0000
0.1429 0.0000 0.2857
>>
1.0000 -0.0000 0.0000
0.0000 1.0000 0.0000
0.0000 0.0000 1.0000
. . .