У меня есть 4 файла .mtx, из которых я читаю значения.Два из них отлично работают при чтении без проблем и выдают правильные результаты в файл .DAT.Тем не менее, последние 2 являются чрезвычайно большими файлами;кажется, что код правильно читает из файлов и работает, но я не получаю никаких выходных данных и никаких ошибок при чтении из этих 2 ... даже таймер кода не печатает время.Любая помощь высоко ценится!Вот код:
program proj2matrixC40
implicit none
integer,parameter::dp=selected_real_kind(15,307)
! Set Global Variables
real(kind=dp), allocatable::Ax(:,:),A(:,:),Iglobal(:,:)
integer::At(1,3)
integer::nnz,w,n,k,ii,ff,kk
real(kind=dp)::t1,t2
call cpu_time(t1)
open(unit=78,file="e40r5000.mtx",status='old')
read(78,*) At
close(unit=78)
nnz = At(1,3)
n = At(1,1)
k = 40
kk = 35
allocate(Ax(nnz+1,3),A(nnz,3),Iglobal(k,k))
open(unit=61,file="e40r5000.mtx",status='old')
do w=1,nnz+1
read(61,*) Ax(w,:)
end do
open (unit = 53, file = "proj2matrixC40points.dat")
do ff=1,k
do ii=1,k
Iglobal(ii,ff) = (ii/ff)*(ff/ii)
end do
end do
A(1:nnz,:) = Ax(2:nnz+1,:)
call Arno(A)
call cpu_time(t2)
print '("Time elapsed = ",f10.8," seconds")', (t2 - t1)
contains
subroutine Arno(a)
real(kind=dp), intent(in)::a(:,:)
real(kind=dp),dimension(k,k)::H
real(kind=dp),dimension(k,k+1)::u,q,qconj
real(kind=dp),dimension(k,1)::x0
integer::j,f
call random_number(x0)
q(:,1) = x0(:,1)/norm2(x0(:,1))
do f=1,k
call spmat(a,q(:,f),u(:,f))
do j=1,f
qconj(j,:) = (q(:,j))
H(j,f) = dot_product(qconj(j,:),u(:,f))
u(:,f) = u(:,f) - H(j,f)*q(:,j)
end do
if (f.lt.k) then
H(f+1,f) = norm2(u(:,f))
if (H(f+1,f)==0) then
print *, "Matrix is reducible"
stop
end if
q(:,f+1) = u(:,f)/H(f+1,f)
end if
if (f==k) then
call qrit(H)
end if
end do
end subroutine
! QR Iteration with Shifts Subroutine
subroutine qrit(a)
real(kind=dp), intent(in)::a(:,:)
real(kind=dp)::sigmak
real(kind=dp),dimension(kk,k)::dia
real(kind=dp),dimension(k,k)::Qfinal,Rfinal,HH
real(kind=dp),dimension(k,k,kk)::H0,needQR
integer::v,z
HH = a
H0(:,:,1) = HH
do v=1,kk
sigmak = H0(k,k,v)
if (v-1==0) then
needQR(:,:,v) = HH - sigmak*Iglobal
else
needQR(:,:,v) = H0(:,:,v-1) - sigmak*Iglobal
end if
call givens2(needQR(:,:,v),Rfinal,Qfinal)
H0(:,:,v) = matmul(Rfinal,Qfinal) + sigmak*Iglobal
do z = 1,k
dia(v,z) = H0(z,z,v)
write(53,*) v," ", dia(v,z) ! Write values to .DAT file
end do
end do
end subroutine
! Sparse Matrix Vector Multiplication Subroutine
subroutine spmat(a,b,c)
real(kind=dp), intent(in)::a(:,:)
real(kind=dp), intent(in), dimension(k,1)::b
real(kind=dp), intent(out), dimension(k,1)::c
integer::m,rowi,columni
real(kind=dp), dimension(k,1)::x,y
x = b
y(:,1) = 0
do m = 1,nnz
rowi = a(m,1)
columni = a(m,2)
y(rowi,1) = y(rowi,1) + a(m,3)*x(columni,1)
end do
c(:,1) = y(:,1)
end subroutine
! QR Factorization Givens Rotations Subroutine
subroutine givens2(a,Rfinal,Qfinal)
real(kind=dp), intent(in)::a(:,:)
real(kind=dp), dimension(k,k,(k*k))::G,QQ
real(kind=dp), dimension(k,k), intent(out)::Rfinal,Qfinal
real(kind=dp), dimension(k,k)::I2,y,aa
real(kind=dp), dimension(1,k)::ek1,ek2
real(kind=dp)::c,s
integer::kt,m,nn,j,i,l,p
m = size(a,1)
nn = size(a,2)
aa = a
i = 1
do kt=1,nn-1
do j=m,kt+1,-1
if (aa(j,kt).eq.0) then
continue
else
ek1(1,:) = 0
ek2(1,:) = 0
do p=1,m
do l=1,m
I2(l,p) = (l/p)*(p/l)
end do
end do
c = aa(kt,kt)/sqrt(aa(kt,kt)**2 + aa(j,kt)**2)
s = aa(j,kt)/sqrt(aa(kt,kt)**2 + aa(j,kt)**2)
ek1(1,kt) = c
ek1(1,j) = s
ek2(1,kt) = -s
ek2(1,j) = c
I2(kt,:) = ek1(1,:)
I2(j,:) = ek2(1,:)
G(:,:,i) = I2
if (i.eq.1) then
QQ(:,:,i) = G(:,:,i)
else
QQ(:,:,i) = matmul(G(:,:,i),QQ(:,:,i-1))
end if
y = matmul(G(:,:,i),aa)
aa = y
if (kt.eq.nn-1) then
if (j.eq.kt+1) then
Qfinal = transpose(QQ(:,:,i))
Rfinal = aa
end if
end if
i = i + 1
end if
end do
end do
end subroutine
end program proj2matrixC40
Пара замечаний.Строка, вокруг которой я ставлю звездочки (для этого вопроса) call mat_print('H',H)
, не может быть удалена, в противном случае я получаю неправильные ответы (это странно ... мысли?).Кроме того, чтобы ваш компьютер не зависал при открытии больших файлов, они называются 'e40r5000.mtx' и 's3dkt3m2.mtx' (это две проблемы, с которыми у меня проблемы).Я использую gfortran версии 8.1.0
Вот ссылка на файлы
https://1drv.ms/f/s!AjG0dE43DVddaJfY62ABE8Yq3CI