Line data Source code
1 0 : subroutine estrad(n,nerror,fjac,ldjac,mljac,mujac,fmas,ldmas,mlmas,mumas,
2 0 : & h,dd1,dd2,dd3,fcn,nfcn,y0,y,ijob,x,m1,m2,nm1,
3 0 : & e1_1D,lde1,z1,z2,z3,cont,f1,f2,ip1,iphes,scal,err,
4 : & first,reject,fac1,lrpar,rpar,lipar,ipar,ier,
5 : & mle,mue,mbjac,mbb,mdiag,mdiff,mbdiag,
6 0 : & decsol,decsols,nzmax,isparse,ia,ja,sa,
7 : & lrd,rpar_decsol,lid,ipar_decsol)
8 : implicit none
9 : interface
10 : #include "num_fcn.dek"
11 : #include "mtx_decsol.dek"
12 : #include "mtx_decsols.dek"
13 : end interface
14 : integer, pointer :: ip1(:) ! (nm1)
15 : integer :: n, iphes(n), nerror, ldjac, mljac, mujac, ldmas, mlmas, mumas
16 : integer :: nfcn, ijob, m1, m2, nm1, lde1, lrpar, lipar, ier, mle, mue
17 : integer :: mbjac, mbb, mdiag, mdiff, mbdiag, nzmax, isparse, lrd, lid
18 : integer :: ia(:) ! (n+1)
19 : integer :: ja(:) ! (nzmax)
20 : double precision :: sa(nzmax)
21 : real(dp), intent(inout), pointer :: rpar_decsol(:) ! (lrd)
22 : integer, intent(inout), pointer :: ipar_decsol(:) ! (lid)
23 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
24 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
25 :
26 : double precision :: fjac(ldjac,n), fmas(ldmas,nm1)
27 : double precision, pointer :: e1_1D(:)
28 : double precision :: x, scal(n), y0(n), y(n)
29 : double precision, pointer, dimension(:) :: cont, z1, z2, z3, f1, f2
30 : double precision :: dd1, dd2, dd3, h, fac1, err
31 : logical :: first, reject
32 :
33 : ! LOCALS
34 0 : double precision :: hee1, hee2, hee3, sum, sum1
35 : integer :: i, mm, j, k, im1
36 : real(dp), pointer :: e1(:,:)
37 : real(dp), pointer, dimension(:) :: p1
38 :
39 0 : e1(1:lde1,1:nm1) => e1_1D(1:lde1*nm1)
40 0 : p1(1:n) => cont(m1+1:m1+n)
41 :
42 0 : ier = 0
43 0 : hee1=dd1/h
44 0 : hee2=dd2/h
45 0 : hee3=dd3/h
46 :
47 0 : mm=m1/m2
48 0 : GOTO (1,2,3,4,5,6,55,8,9,55,11,12,13,14,15), ijob
49 : !
50 : 1 continue
51 : ! ------ b=identity, jacobian a full matrix
52 0 : do i=1,n
53 0 : f2(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
54 0 : cont(i)=f2(i)+y0(i)
55 : end do
56 0 : call decsol(1,n,lde1,e1_1D,n,n,cont,ip1,lrd,rpar_decsol,lid,ipar_decsol,ier)
57 0 : GOTO 77
58 : !
59 : 8 continue
60 : ! ------ b=identity, jacobian a sparse matrix
61 0 : do i=1,n
62 0 : f2(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
63 0 : cont(i)=f2(i)+y0(i)
64 : end do
65 0 : call decsols(1,n,nzmax,ia,ja,sa,cont,lrd,rpar_decsol,lid,ipar_decsol,ier)
66 0 : GOTO 77
67 : !
68 : 11 continue
69 : ! ------ b=identity, jacobian a full matrix, second order
70 0 : do i=1,n
71 0 : f2(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
72 0 : cont(i)=f2(i)+y0(i)
73 : end do
74 : 48 continue
75 0 : do j=1,m2
76 0 : sum1=0.d0
77 0 : do k=mm-1,0,-1
78 0 : sum1=(cont(j+k*m2)+sum1)/fac1
79 0 : do i=1,nm1
80 0 : im1=i+m1
81 0 : cont(im1)=cont(im1)+fjac(i,j+k*m2)*sum1
82 : end do
83 : end do
84 : end do
85 0 : call decsol(1,nm1,lde1,e1_1D,nm1,nm1,p1,ip1,lrd,rpar_decsol,lid,ipar_decsol,ier)
86 0 : do i=m1,1,-1
87 0 : cont(i)=(cont(i)+cont(m2+i))/fac1
88 : end do
89 : GOTO 77
90 : !
91 : 2 continue
92 : ! ------ b=identity, jacobian a banded matrix
93 0 : do i=1,n
94 0 : f2(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
95 0 : cont(i)=f2(i)+y0(i)
96 : end do
97 0 : call decsol(1,n,lde1,e1_1D,mle,mue,cont,ip1,lrd,rpar_decsol,lid,ipar_decsol,ier)
98 0 : GOTO 77
99 : !
100 : 12 continue
101 : ! ------ b=identity, jacobian a banded matrix, second order
102 0 : do i=1,n
103 0 : f2(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
104 0 : cont(i)=f2(i)+y0(i)
105 : end do
106 : 45 continue
107 0 : do j=1,m2
108 0 : sum1=0.d0
109 0 : do k=mm-1,0,-1
110 0 : sum1=(cont(j+k*m2)+sum1)/fac1
111 0 : do i=max(1,j-mujac),min(nm1,j+mljac)
112 0 : im1=i+m1
113 0 : cont(im1)=cont(im1)+fjac(i+mujac+1-j,j+k*m2)*sum1
114 : end do
115 : end do
116 : end do
117 0 : call decsol(1,nm1,lde1,e1_1D,mle,mue,p1,ip1,lrd,rpar_decsol,lid,ipar_decsol,ier)
118 0 : do i=m1,1,-1
119 0 : cont(i)=(cont(i)+cont(m2+i))/fac1
120 : end do
121 : GOTO 77
122 : !
123 : 3 continue
124 : ! ------ b is a banded matrix, jacobian a full matrix
125 0 : do i=1,n
126 0 : f1(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
127 : end do
128 0 : do i=1,n
129 0 : sum=0.d0
130 0 : do j=max(1,i-mlmas),min(n,i+mumas)
131 0 : sum=sum+fmas(i-j+mbdiag,j)*f1(j)
132 : end do
133 0 : f2(i)=sum
134 0 : cont(i)=sum+y0(i)
135 : end do
136 0 : call decsol(1,n,lde1,e1_1D,n,n,cont,ip1,lrd,rpar_decsol,lid,ipar_decsol,ier)
137 0 : GOTO 77
138 : !
139 : 9 continue
140 : ! ------ b is a banded matrix, jacobian a sparse matrix
141 0 : do i=1,n
142 0 : f1(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
143 : end do
144 0 : do i=1,n
145 0 : sum=0.d0
146 0 : do j=max(1,i-mlmas),min(n,i+mumas)
147 0 : sum=sum+fmas(i-j+mbdiag,j)*f1(j)
148 : end do
149 0 : f2(i)=sum
150 0 : cont(i)=sum+y0(i)
151 : end do
152 0 : call decsols(1,n,nzmax,ia,ja,sa,z1,lrd,rpar_decsol,lid,ipar_decsol,ier)
153 0 : GOTO 77
154 : !
155 : 13 continue
156 : ! ------ b is a banded matrix, jacobian a full matrix, second order
157 0 : do i=1,m1
158 0 : f2(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
159 0 : cont(i)=f2(i)+y0(i)
160 : end do
161 0 : do i=m1+1,n
162 0 : f1(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
163 : end do
164 0 : do i=1,nm1
165 0 : sum=0.d0
166 0 : do j=max(1,i-mlmas),min(nm1,i+mumas)
167 0 : sum=sum+fmas(i-j+mbdiag,j)*f1(j+m1)
168 : end do
169 0 : im1=i+m1
170 0 : f2(im1)=sum
171 0 : cont(im1)=sum+y0(im1)
172 : end do
173 : GOTO 48
174 : !
175 : 4 continue
176 : ! ------ b is a banded matrix, jacobian a banded matrix
177 0 : do i=1,n
178 0 : f1(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
179 : end do
180 0 : do i=1,n
181 0 : sum=0.d0
182 0 : do j=max(1,i-mlmas),min(n,i+mumas)
183 0 : sum=sum+fmas(i-j+mbdiag,j)*f1(j)
184 : end do
185 0 : f2(i)=sum
186 0 : cont(i)=sum+y0(i)
187 : end do
188 0 : call decsol(1,n,lde1,e1_1D,mle,mue,cont,ip1,lrd,rpar_decsol,lid,ipar_decsol,ier)
189 0 : GOTO 77
190 : !
191 : 14 continue
192 : ! ------ b is a banded matrix, jacobian a banded matrix, second order
193 0 : do i=1,m1
194 0 : f2(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
195 0 : cont(i)=f2(i)+y0(i)
196 : end do
197 0 : do i=m1+1,n
198 0 : f1(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
199 : end do
200 0 : do i=1,nm1
201 0 : sum=0.d0
202 0 : do j=max(1,i-mlmas),min(nm1,i+mumas)
203 0 : sum=sum+fmas(i-j+mbdiag,j)*f1(j+m1)
204 : end do
205 0 : im1=i+m1
206 0 : f2(im1)=sum
207 0 : cont(im1)=sum+y0(im1)
208 : end do
209 : GOTO 45
210 : !
211 : 5 continue
212 : ! ------ b is a full matrix, jacobian a full matrix
213 0 : do i=1,n
214 0 : f1(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
215 : end do
216 0 : do i=1,n
217 : sum=0.d0
218 0 : do j=1,n
219 0 : sum=sum+fmas(i,j)*f1(j)
220 : end do
221 0 : f2(i)=sum
222 0 : cont(i)=sum+y0(i)
223 : end do
224 0 : call decsol(1,n,lde1,e1_1D,n,n,cont,ip1,lrd,rpar_decsol,lid,ipar_decsol,ier)
225 0 : GOTO 77
226 : !
227 : 15 continue
228 : ! ------ b is a banded matrix, jacobian a full matrix, second order
229 0 : do i=1,m1
230 0 : f2(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
231 0 : cont(i)=f2(i)+y0(i)
232 : end do
233 0 : do i=m1+1,n
234 0 : f1(i)=hee1*z1(i)+hee2*z2(i)+hee3*z3(i)
235 : end do
236 0 : do i=1,nm1
237 : sum=0.d0
238 0 : do j=1,nm1
239 0 : sum=sum+fmas(i,j)*f1(j+m1)
240 : end do
241 0 : im1=i+m1
242 0 : f2(im1)=sum
243 0 : cont(im1)=sum+y0(im1)
244 : end do
245 : GOTO 48
246 : !
247 : 6 continue
248 : ! ------ b is a full matrix, jacobian a banded matrix
249 : ! ------ this option is not provided
250 0 : return
251 : !
252 : ! --------------------------------------
253 : !
254 : 77 continue
255 0 : err=0.d0
256 0 : do i=1,n
257 0 : err=err+(cont(i)/scal(i))**2
258 : end do
259 0 : err=max(sqrt(err/n),1.d-10)
260 : !
261 0 : if (err<1.d0) return
262 0 : if (first.or.reject) then
263 0 : do i=1,n
264 0 : cont(i)=y(i)+cont(i)
265 : end do
266 0 : call fcn(n,x,h,cont,f1,lrpar,rpar,lipar,ipar,ier)
267 0 : if (ier /= 0) return
268 0 : nfcn=nfcn+1
269 0 : do i=1,n
270 0 : cont(i)=f1(i)+f2(i)
271 : end do
272 0 : GOTO (31,32,31,32,31,32,55,55,33,55,41,42,41,42,41), ijob
273 : ! ------ full matrix option
274 : 31 continue
275 0 : call decsol(1,n,lde1,e1_1D,n,n,cont,ip1,lrd,rpar_decsol,lid,ipar_decsol,ier)
276 0 : GOTO 88
277 : ! ------ full matrix option, second order
278 : 41 continue
279 0 : do j=1,m2
280 0 : sum1=0.d0
281 0 : do k=mm-1,0,-1
282 0 : sum1=(cont(j+k*m2)+sum1)/fac1
283 0 : do i=1,nm1
284 0 : im1=i+m1
285 0 : cont(im1)=cont(im1)+fjac(i,j+k*m2)*sum1
286 : end do
287 : end do
288 : end do
289 0 : call decsol(1,nm1,lde1,e1_1D,nm1,nm1,p1,ip1,lrd,rpar_decsol,lid,ipar_decsol,ier)
290 0 : do i=m1,1,-1
291 0 : cont(i)=(cont(i)+cont(m2+i))/fac1
292 : end do
293 : GOTO 88
294 : ! ------ banded matrix option
295 : 32 continue
296 0 : call decsol(1,n,lde1,e1_1D,mle,mue,cont,ip1,lrd,rpar_decsol,lid,ipar_decsol,ier)
297 0 : GOTO 88
298 : ! ------ banded matrix option, second order
299 : 42 continue
300 0 : do j=1,m2
301 0 : sum1=0.d0
302 0 : do k=mm-1,0,-1
303 0 : sum1=(cont(j+k*m2)+sum1)/fac1
304 0 : do i=max(1,j-mujac),min(nm1,j+mljac)
305 0 : im1=i+m1
306 0 : cont(im1)=cont(im1)+fjac(i+mujac+1-j,j+k*m2)*sum1
307 : end do
308 : end do
309 : end do
310 0 : call decsol(1,nm1,lde1,e1_1D,mle,mue,p1,ip1,lrd,rpar_decsol,lid,ipar_decsol,ier)
311 0 : do i=m1,1,-1
312 0 : cont(i)=(cont(i)+cont(m2+i))/fac1
313 : end do
314 : GOTO 88
315 : ! ------ sparse matrix option
316 : 33 continue
317 0 : call decsols(1,n,nzmax,ia,ja,sa,cont,lrd,rpar_decsol,lid,ipar_decsol,ier)
318 0 : GOTO 88
319 : ! -----------------------------------
320 : 88 continue
321 0 : err=0.d0
322 0 : do i=1,nerror
323 0 : err=err+(cont(i)/scal(i))**2
324 : end do
325 0 : err=max(sqrt(err/nerror),1.d-10)
326 : end if
327 : return
328 : ! -----------------------------------------------------------
329 : 55 continue
330 0 : write(*,*) 'estrad: invalid ijob', ijob
331 0 : call mesa_error(__FILE__,__LINE__) ! estrad
332 : end subroutine estrad
|
