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
|
program phitest
c-------------------------------------------------------------------
c Test program for exponential propagator using Arnoldi approach
c This main program is a very simple test using diagonal matrices
c (Krylov subspace methods are blind to the structure of the matrix
c except for symmetry). This provides a good way of testing the
c accuracy of the method as well as the error estimates. This test
c program tests the phi-function variant instead of the exponential.
c
c The subroutine phipro which is tested here computes an approximation
c to the vector
c
c w(tn) = w(t0) + tn * phi( - A * tn ) * (r - A w(t0))
c
c where phi(z) = (1-exp(z)) / z
c
c i.e. it solves dw/dt = -A w + r in [t0,t0+ tn] (returns only w(t0+tn))
c
c In this program t0=0, tn is input by the user and A is a simple
c diagonal matrix.
c
c-------------------------------------------------------------------
implicit real*8 (a-h,o-z)
parameter (nmax = 400, ih0=60, ndmx=20,nzmax = 7*nmax)
real*8 a(nzmax), u(ih0*nmax), w(nmax),w1(nmax),x(nmax),y(nmax),
* r(nmax)
integer ioff(10)
data iout/6/,a0/0.0/,b0/1.0/,epsmac/1.d-10/,eps/1.d-10/
c
c dimendion of matrix
c
n = 100
c---------------------------define matrix -----------------------------
c A is a single diagonal matrix (ndiag = 1 and ioff(1) = 0 )
c-----------------------------------------------------------------------
ndiag = 1
ioff(1) = 0
c
c entriesin the diagonal are uniformly distributed.
c
h = 1.0d0 / real(n+1)
do 1 j=1, n
a(j) = real(j+1)*h
1 continue
c--------
write (6,'(10hEnter tn: ,$)')
read (5,*) tn
c
write (6,'(36hEpsilon (desired relative accuracy): ,$)')
read (5,*) eps
c-------
write (6,'(36h m (= dimension of Krylov subspace): ,$)')
read (5,*) m
c
c define initial conditions to be (1,1,1,1..1)^T
c
do 2 j=1,n
w(j) = 1.0 - 1.0/real(j+10)
r(j) = 1.0
w1(j) = w(j)
2 continue
c
c
c HERE IT IS DONE IN 10 SUBSTEPS
c nsteps = 10
c tnh = tn/real(nsteps)
c MARCHING LOOP
c do jj =1, nsteps
c call phiprod (n, m, eps, tnh, u, w, r, x, y, a, ioff, ndiag)
call phiprod (n, m, eps, tn, u, w, r, x, y, a, ioff, ndiag)
c enddo
print *, ' final answer '
print *, (w(k),k=1,20)
c
c exact solution
c
do 4 k=1,n
w1(k)=w1(k)+((1.0 - dexp(-a(k)*tn) )/a(k))*(r(k)
+ - a(k)*w1(k) )
4 continue
c
c exact solution
c
print *, ' exact answer '
print *, (w1(k),k=1,20)
c
c computing actual 2-norm of error
c
t = 0.0d0
do 47 k=1,n
t = t+ (w1(k)-w(k))**2
47 continue
t = dsqrt(t / ddot(n, w,1,w,1) )
c
write (6,*) ' final error', t
c-----------------------------------------------------------------------
stop
end
c-----------------------------------------------------------------------
subroutine oped(n,x,y,diag,ioff,ndiag)
c======================================================
c this kernel performs a matrix by vector multiplication
c for a diagonally structured matrix stored in diagonal format
c======================================================
implicit real*8 (a-h,o-z)
real*8 x(n), y(n), diag(n,ndiag)
common nope, nmvec
integer j, n, ioff(ndiag)
CDIR$ IVDEP
do 1 j=1, n
y(j) = 0.00
1 continue
c
do 10 j=1,ndiag
io = ioff(j)
i1=max0(1,1-io)
i2=min0(n,n-io)
CDIR$ IVDEP
do 9 k=i1,i2
y(k) = y(k)+diag(k,j)*x(k+io)
9 continue
10 continue
nmvec = nmvec + 1
nope = nope + 2*ndiag*n
return
end
c
|
