1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
subroutine worker
use globals
use auxiliary
implicit none
logical :: conv
integer :: seed,ierr,recvd_tag,tag
integer, dimension(mpi_status_size) :: status
double precision :: eig
double precision, allocatable, dimension(:,:) :: A
! Receive matrix size
call mpi_bcast(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
allocate(A(n,n))
do while(.true.)
call mpi_recv(seed,1,MPI_INTEGER,0,MPI_ANY_TAG,MPI_COMM_WORLD,STATUS,ierr)
recvd_tag = status (MPI_TAG)
if(recvd_tag.eq.EXIT_TAG) then ! The exit tag was received, wind down and exit.
print *,'Worker ',proc_num,' exiting...'
exit
end if
! The seed was received and more data are requested, so compute an eigenvalue:
call open_random_stream(seed)
call make_random_matrix(A,seed)
call compute_eig(A,eig,conv)
! Send result
if(conv) then
tag=0
else
tag=1
end if
call mpi_send(eig,1,MPI_DOUBLE_PRECISION,0,tag,MPI_COMM_WORLD,ierr)
call close_random_stream()
end do
deallocate(A)
return
end subroutine worker
subroutine make_random_matrix(A,seed)
use globals
use auxiliary
implicit none
integer :: seed,ierr,i,j
double precision :: A(n,n)
do i=1,n
ierr=vdRNGGaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF,stream,1,A(i,i),0d0,2d0)
do j=1,i-1
ierr=vdRNGGaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF,stream,1,A(i,j),0d0,1d0)
A(j,i)=A(i,j)
end do
end do
return
end subroutine make_random_matrix
subroutine compute_eig(A,eig,conv)
use globals
use auxiliary
implicit none
logical :: conv
integer :: ierr,M,lwork,info,i
double precision :: A(n,n),eig,eps,mu,wr(n),wi(n),B(n,n),vl(1,n),vr(1,n),eig2
double precision, allocatable, dimension(:) :: work
lwork=4*n
allocate(work(lwork)) ! Bit of memory needed by LAPACK
eps=1d-5 ! Relative convergence criterion for power iteration (bad practice to hard-code this!)
M=int(n/2) ! Maximal number of iterations
mu=0d0
call power(A,mu,eps,M,eig,conv)
! Check if the eigenvalue is positive (they are symmetrically distributed about 0). If not, re-compute with shift:
if(conv) then
if(eig.lt.0) then
mu=-eig
call power(A,mu,eps,M,eig,conv)
end if
end if
! If power iteration fails to converge, compute the whole spectrum with Shur decomposition:
if(.not.(conv)) then
call DGEEV('N','N',n,A,n,wr,wi,vl,1,vr,1,work,lwork,info)
if(info.eq.0) then
eig=wr(maxloc(wr,1)) ! Select the largest eigenvalue
conv=.true.
else
print *,'error in DGEEV ',info
end if
end if
deallocate(work)
return
end subroutine compute_eig
subroutine power(A,mu,eps,M,eig,conv)
use globals
use auxiliary
implicit none
logical :: conv
integer :: i,M,ierr,j
double precision :: A(n,n),v(n),y(n),eig,eps,mu,err,nrm,ndp,mem
conv=.false.
! Make random initial vector
ierr=vdrnggaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF,stream,n,v,0d0,1d0)
v=v/sqrt(sum(v**2))
mem=-100d0
do i=1,M
call dgemv('N',n,n,1d0,A,n,v,1,0d0,y,1)
y=y+mu*v
nrm=sqrt(sum(y**2))
eig=dot_product(v,y) ! Rayleigh quotient approximates the eigenvalue
eig=eig-mu
err=abs(mem-eig)/abs(mem)
if(err.lt.eps) then
conv=.true.
exit
end if
v=y/nrm
mem=eig
end do
return
end subroutine power
|
