Line data Source code
1 80 : subroutine vintk(x,a,y,ii,nn,idima,idimy,karg)
2 : c
3 : c 2k-1 -th order integration routine
4 : c ***********************************
5 : c
6 : c sets integral of a(i,n) with respect to x into y(i,n),i=1,ii,n=1, nn.
7 : c the values of y(i,1), i = 1, ii, must be supplied.
8 : c idima and idimy are first dimensions of a and y respectively.
9 : c if x is found to be non-monotonic ii is set to 0.
10 : c
11 : c if vintk is called with karg .le. 0, it is assumed that a table
12 : c of integration coefficients has been set up.
13 : c
14 : c vintk uses as work space
15 : c wrk(2k*2k)
16 : c the size of wrk set in vintk is 100, which is sufficient, say, for
17 : c 9-th order integration in 20 dependent variables. if more
18 : c work space is needed it must be set in common/wrklir/
19 : c
20 : c note that this version of vintk also uses storage of size
21 : c
22 : c dc(2*k*(nn-1))
23 : c
24 : c in common /cintst/ to store integration coefficients. the
25 : c default size of 500 probably in general has to be increased in
26 : c the calling programme.
27 : c
28 : c original version: 19/3/1986
29 : c
30 : c Double precision version.
31 : c ++++++++++++++++++++++++
32 : c
33 : c Dated 10/3/90.
34 : c
35 : c Modified 26/4/99, fixing redefinition of k when nn is too small
36 : c
37 : c Note: this double precision version of the routine has same name
38 : c as single precision version, but is distinguished by its file name
39 : c
40 : implicit double precision(a-h,o-z)
41 : dimension x(1),a(idima,1),y(idimy,1)
42 : common/wrklir/ wrk(100)
43 : common/cintst/ dc(500)
44 : c
45 : c common defining standard input and output
46 : c
47 : common/cstdio/ istdin, istdou, istdpr, istder
48 : c
49 : save /cintst/,k,k2
50 : c
51 : c test for setting up coefficients
52 : c
53 80 : if(karg.le.0) go to 60
54 : c
55 1 : k=karg
56 : c
57 : c check order
58 : c
59 1 : k2=2*k
60 1 : if(nn.lt.k2) then
61 0 : if(istdpr.gt.0) write(istdpr,100) k2,nn
62 0 : k=nn/2
63 0 : k2=2*k
64 : end if
65 : c
66 : c first dimension of wrk, storage information
67 : c
68 1 : iw=k2
69 1 : mst=iw*(iw-1)+2
70 : c
71 : c total range
72 : c
73 1 : diff=x(nn)-x(1)
74 1 : if(diff.eq.0) go to 95
75 : c
76 : c start loop for setting coefficients
77 : c
78 1 : 20 jdc=0
79 : c
80 1468 : do 50 n=2,nn
81 1467 : n1=n-1
82 1467 : xtr=x(n1)
83 : c
84 : c determine i
85 : c
86 1467 : i=min0(max0(0,n1-k),nn-k2)
87 : c length of interval
88 1467 : dlx=x(i+k2)-x(i+1)
89 1467 : if(diff*dlx.le.0) go to 97
90 : c
91 : c set equations
92 : c
93 1467 : 32 dlxi=1.d0/dlx
94 : c
95 1467 : m=0
96 7335 : do 35 l=1,k2
97 5868 : l1=i+l
98 5868 : xl=(x(l1)-xtr)*dlxi
99 5868 : aa=1.d0
100 5868 : m=m+1
101 5868 : wrk(m)=aa
102 24939 : do 35 ir=2,k2
103 17604 : m=m+1
104 17604 : aa=xl*aa
105 23472 : 35 wrk(m)=aa
106 : c
107 : c set right hand sides
108 : c
109 1467 : j1=jdc
110 1467 : xl=(x(n)-xtr)*dlxi
111 1467 : aa=1
112 : c
113 7335 : do 40 j=1,k2
114 5868 : aa=aa*xl
115 5868 : j1=j1+1
116 7335 : 40 dc(j1)=dlx*aa/j
117 : c
118 : c solve equations
119 : c
120 1467 : call leq(wrk,dc(jdc+1),k2,1,iw,1,det)
121 : c
122 1467 : if(det.eq.0) go to 98
123 : c
124 1468 : 50 jdc=jdc+k2
125 : c
126 : c end setting coefficients
127 : c
128 1 : if(istdpr.gt.0) write(istdpr,110) jdc
129 : c
130 : c *****************************************
131 : c
132 : c set integrals
133 : c
134 80 : 60 jdc=0
135 : c
136 117440 : do 70 n=2,nn
137 117360 : n1=n-1
138 : c
139 : c determine i
140 : c
141 117360 : i=min0(max0(0,n1-k),nn-k2)
142 : c
143 369684 : do 65 is=1,ii
144 252324 : j1=jdc
145 252324 : sum=0
146 1261620 : do 62 l=1,k2
147 1009296 : j1=j1+1
148 1261620 : 62 sum=sum+dc(j1)*a(is,i+l)
149 : c
150 369684 : 65 y(is,n)=y(is,n1)+sum
151 : c
152 117440 : 70 jdc=jdc+k2
153 : c
154 80 : return
155 : c
156 : c diagnostics
157 : c
158 0 : 95 if(istdpr.gt.0) write(istdpr,120) x(1)
159 0 : ii=0
160 0 : return
161 0 : 97 i1=i+1
162 0 : ik2=i+k2
163 0 : if(istdpr.gt.0) write(istdpr,130) x(1),x(nn),i1,x(i1),ik2,x(ik2)
164 0 : ii=0
165 0 : return
166 0 : 98 i1=i+1
167 0 : ik2=i+k2
168 0 : if(istdpr.gt.0) write(istdpr,140) (l,x(l),l=i1,ik2)
169 0 : ii=0
170 0 : return
171 : 100 format(1x,10('*'),' 2k+1 =',i4,' is greater than nn =',i4,
172 : * ' in vintk.'/11x,'k has been reset to nn/2')
173 : 110 format(//' storage needed in common/cintst/ in vintk:',
174 : * i7)
175 : 120 format(//1x,10('*'),' range is zero in vintk. x(1)=',
176 : * 1pe15.5,' = x(nn)'//)
177 : 130 format(//1x,10('*'),' independent variable is not monotonic',
178 : * ' in vintk'//' x(1)=',1pe15.5,' x(nn)=',e15.5,
179 : * 2(' x(',i4,')=',e15.5)//)
180 : 140 format(//1x,10('*'),' points coincide in vintk'/
181 : * (5(' x(',i4,')=',1pe13.4)))
182 : end
|
