Line data Source code
1 162344 : subroutine slvrod(n,fjac,ldjac,mljac,mujac,fmas,ldmas,mlmas,mumas,
2 81172 : & m1,m2,nm1,fac1,e_1D,lde,ip,dy,ak,fx,ynew,hd,ijob,not_stage1,
3 : & mle,mue,mbjac,mbb,mdiag,mdiff,mbdiag,
4 : & decsol,decsols,decsolblk,
5 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk1, uf_dblk1, uf_ublk1,
6 81172 : & nzmax,isparse,ia,ja,sa,lrd,rpar_decsol,lid,ipar_decsol,ier)
7 : use mtx_lib
8 : implicit none
9 : interface
10 : #include "mtx_decsol.dek"
11 : #include "mtx_decsols.dek"
12 : #include "mtx_decsolblk.dek"
13 : end interface
14 : integer, intent(in) :: caller_id, nvar, nz
15 : real(dp), dimension(:), pointer, intent(inout) :: lblk, dblk, ublk ! =(nvar,nvar,nz)
16 : real(dp), pointer, dimension(:) :: uf_lblk1, uf_dblk1, uf_ublk1
17 : integer :: ia(:) ! (n+1)
18 : integer :: ja(:) ! (nzmax)
19 : real(dp) :: sa(:) ! (nzmax)
20 : integer :: mle, mue, mbjac, mbb, mdiag, mdiff, mbdiag, nzmax, isparse, lrd, lid
21 : integer :: n, ldjac, mljac, mujac, ldmas, mlmas, mumas, m1, m2, nm1, lde, ijob
22 : integer :: ier
23 : real(dp) :: fjac(ldjac,n), fmas(ldmas,nm1), dy(n)
24 : real(dp) :: fac1, hd, fx(n), ynew(n)
25 : logical :: not_stage1
26 : real(dp), intent(inout), pointer :: rpar_decsol(:) ! (lrd)
27 : integer, intent(inout), pointer :: ipar_decsol(:) ! (lid)
28 : integer, pointer :: ip(:) ! (nm1)
29 : real(dp), pointer :: ak(:)
30 : real(dp), pointer :: e_1D(:) ! =(lde,nm1)
31 :
32 : ! LOCALS
33 : integer :: i, j
34 81172 : real(dp) :: sum
35 81172 : real(dp), pointer :: r1(:)
36 162344 : real(dp), target, dimension(nvar*nz) :: x0_ary, b_ary, r_ary, Ax0_ary
37 81172 : real(dp), pointer, dimension(:,:) :: ak2, x0, b, r, Ax0
38 81172 : real(dp), pointer, dimension(:,:,:) :: uf_lblk, uf_dblk, uf_ublk
39 :
40 : !
41 81172 : if (hd == 0.d0) then
42 9513100 : do i=1,n
43 9513100 : ak(i)=dy(i)
44 : end do
45 : else
46 0 : do i=1,n
47 0 : ak(i)=dy(i)+hd*fx(i)
48 : end do
49 : end if
50 :
51 : ! write(*,*) 'slvrod ijob', ijob
52 :
53 : !
54 81172 : if (nvar > 0) then
55 :
56 0 : if (nvar*nz /= n) stop 'bad nvar*nz /= n'
57 :
58 0 : if (not_stage1) then
59 : !do i=1,n
60 : ! write(*,*) 'in ak ynew', i, ak(i), ynew(i)
61 : !end do
62 0 : do i=1,n
63 0 : ak(i)=ak(i)+ynew(i)
64 : end do
65 : end if
66 :
67 : !do i=1,n
68 : ! write(*,*) 'in', i, ak(i)
69 : !end do
70 :
71 0 : ak2(1:nvar,1:nz) => ak(1:n)
72 0 : x0(1:nvar,1:nz) => x0_ary(1:n)
73 0 : b(1:nvar,1:nz) => b_ary(1:n)
74 0 : r(1:nvar,1:nz) => r_ary(1:n)
75 0 : Ax0(1:nvar,1:nz) => Ax0_ary(1:n)
76 :
77 : !do i=1,n
78 : ! b_ary(i) = ak(i) ! save initial rhs
79 : !end do
80 :
81 : ! solve A*ak = b
82 : call decsolblk(
83 0 : & 1,caller_id,nvar,nz,lblk,dblk,ublk,ak,ip,lrd,rpar_decsol,lid,ipar_decsol,ier)
84 : if (ier /= 0) return
85 :
86 :
87 : !do i=1,n
88 : ! write(*,*) 'out', i, ak(i)
89 : !end do
90 : !write(*,*)
91 :
92 0 : return
93 :
94 : do i=1,n
95 : x0_ary(i) = ak(i) ! save initial solution in x0
96 : end do
97 :
98 : ! find dx s.t. A*(x0+dx) = b
99 : ! solve A*dx = b - A*x0
100 : ! refined soln x = x0 + dx
101 :
102 : uf_dblk(1:nvar,1:nvar,1:nz) => uf_dblk1(1:nvar*nvar*nz)
103 : uf_ublk(1:nvar,1:nvar,1:nz) => uf_ublk1(1:nvar*nvar*nz)
104 : uf_lblk(1:nvar,1:nvar,1:nz) => uf_lblk1(1:nvar*nvar*nz)
105 :
106 : call block_dble_mv(nvar,nz,uf_lblk,uf_dblk,uf_ublk,x0,Ax0) ! Ax0 = A*x0
107 :
108 : sum = 0
109 : do i=1,n
110 : r_ary(i) = b_ary(i) - Ax0_ary(i) ! residual
111 : sum = sum + abs(r_ary(i))
112 : end do
113 :
114 : do i=1,n
115 : write(*,*) 'res', i, r_ary(i)
116 : end do
117 : write(*,*)
118 : return
119 :
120 : ! solve A*dx = r
121 : r1(1:n) => r_ary(1:n)
122 : call decsolblk(
123 : & 1,caller_id,nvar,nz,lblk,dblk,ublk,r1,ip,lrd,rpar_decsol,lid,ipar_decsol,ier)
124 : if (ier /= 0) return
125 :
126 : do i=1,n
127 : ak(i) = r_ary(i) + x0_ary(i) ! update solution
128 : end do
129 :
130 : ! check refined solution
131 : write(*,*) 'initial sum', sum
132 :
133 : call block_dble_mv(nvar,nz,uf_lblk,uf_dblk,uf_ublk,ak2,Ax0) ! Ax0 = A*ak
134 :
135 : sum = 0
136 : do i=1,n
137 : r_ary(i) = b_ary(i) - Ax0_ary(i) ! error in solution
138 : sum = sum + abs(r_ary(i))
139 : end do
140 : write(*,*) 'refined sum', sum
141 : write(*,*)
142 :
143 : return
144 : end if
145 :
146 81172 : GOTO (1,2,3,4,5,6,55,8,9,55,11,12,13,13,15), ijob
147 : !
148 : ! -----------------------------------------------------------
149 : !
150 : 1 continue
151 : ! --- b=identity, jacobian a full matrix
152 50988 : if (not_stage1) then
153 185648 : do i=1,n
154 185648 : ak(i)=ak(i)+ynew(i)
155 : end do
156 : end if
157 :
158 : !write(*,*) 'rhs(1)', ak(1)
159 : !write(*,*) 'rhs(2)', ak(2)
160 50988 : call decsol(1,n,lde,e_1D,n,n,ak,ip,lrd,rpar_decsol,lid,ipar_decsol,ier)
161 : !write(*,*) 'ak(1)', ak(1)
162 : !write(*,*) 'ak(2)', ak(2)
163 : !write(*,*)
164 :
165 50988 : return
166 : !
167 : ! -----------------------------------------------------------
168 : !
169 : 8 continue
170 : ! --- b=identity, jacobian a sparse matrix
171 0 : if (not_stage1) then
172 0 : do i=1,n
173 0 : ak(i)=ak(i)+ynew(i)
174 : end do
175 : end if
176 : !write(*,*) 'slvrod call decsols'
177 0 : call decsols(1,n,nzmax,ia,ja,sa,ak,lrd,rpar_decsol,lid,ipar_decsol,ier)
178 : !write(*,*) 'back decsols'
179 0 : return
180 : !
181 : ! -----------------------------------------------------------
182 : !
183 : 9 continue
184 : ! --- b is a banded matrix, jacobian a sparse matrix
185 0 : if (not_stage1) then
186 0 : do i=1,n
187 0 : sum=0.d0
188 0 : do j=max(1,i-mlmas),min(n,i+mumas)
189 0 : sum=sum+fmas(i-j+mbdiag,j)*ynew(j)
190 : end do
191 0 : ak(i)=ak(i)+sum
192 : end do
193 : end if
194 0 : call decsols(1,n,nzmax,ia,ja,sa,ak,lrd,rpar_decsol,lid,ipar_decsol,ier)
195 0 : return
196 : !
197 : ! -----------------------------------------------------------
198 : !
199 : 11 continue
200 : ! --- b=identity, jacobian a full matrix, second order
201 0 : write(*,*) 'SLVROD: second order not supported'
202 0 : call mesa_error(__FILE__,__LINE__)
203 : ! if (not_stage1) then
204 : ! do i=1,n
205 : ! ak(i)=ak(i)+ynew(i)
206 : ! end do
207 : ! end if
208 : ! 48 mm=m1/m2
209 : ! do j=1,m2
210 : ! sum=0.d0
211 : ! do k=mm-1,0,-1
212 : ! jkm=j+k*m2
213 : ! sum=(ak(jkm)+sum)/fac1
214 : ! do i=1,nm1
215 : ! im1=i+m1
216 : ! ak(im1)=ak(im1)+fjac(i,jkm)*sum
217 : ! end do
218 : ! end do
219 : ! end do
220 : ! ptr(1:n) => ak(m1+1:m1+n)
221 : ! call decsol(1,nm1,lde,e_1D,nm1,nm1,ptr,ip,lrd,rpar_decsol,lid,ipar_decsol,ier)
222 : ! do i=m1,1,-1
223 : ! ak(i)=(ak(i)+ak(m2+i))/fac1
224 : ! end do
225 : ! return
226 : !
227 : ! -----------------------------------------------------------
228 : !
229 : 2 continue
230 : ! --- b=identity, jacobian a banded matrix
231 30184 : if (not_stage1) then
232 : !do i=1,n
233 : ! write(*,*) 'in ak ynew', i, ak(i), ynew(i)
234 : !end do
235 7161033 : do i=1,n
236 7161033 : ak(i)=ak(i)+ynew(i)
237 : end do
238 : end if
239 :
240 : !do i=1,n
241 : ! write(*,*) 'in', i, ak(i)
242 : !end do
243 :
244 30184 : call decsol(1,n,lde,e_1D,mle,mue,ak,ip,lrd,rpar_decsol,lid,ipar_decsol,ier)
245 :
246 : !do i=1,n
247 : ! write(*,*) 'out', i, ak(i)
248 : !end do
249 : !write(*,*)
250 30184 : return
251 : !
252 : ! -----------------------------------------------------------
253 : !
254 : 12 continue
255 : ! --- b=identity, jacobian a banded matrix, second order
256 0 : write(*,*) 'SLVROD: second order not supported'
257 0 : call mesa_error(__FILE__,__LINE__)
258 : ! if (not_stage1) then
259 : ! do i=1,n
260 : ! ak(i)=ak(i)+ynew(i)
261 : ! end do
262 : ! end if
263 : ! 45 mm=m1/m2
264 : ! do j=1,m2
265 : ! sum=0.d0
266 : ! do k=mm-1,0,-1
267 : ! jkm=j+k*m2
268 : ! sum=(ak(jkm)+sum)/fac1
269 : ! do i=max(1,j-mujac),min(nm1,j+mljac)
270 : ! im1=i+m1
271 : ! ak(im1)=ak(im1)+fjac(i+mujac+1-j,jkm)*sum
272 : ! end do
273 : ! end do
274 : ! end do
275 : ! ptr(1:n) => ak(m1+1:m1+n)
276 : ! call decsol(1,nm1,lde,e_1D,mle,mue,ptr,ip,lrd,rpar_decsol,lid,ipar_decsol,ier)
277 : ! do i=m1,1,-1
278 : ! ak(i)=(ak(i)+ak(m2+i))/fac1
279 : ! end do
280 : ! return
281 : !
282 : ! -----------------------------------------------------------
283 : !
284 : 3 continue
285 : ! --- b is a banded matrix, jacobian a full matrix
286 0 : if (not_stage1) then
287 0 : do i=1,n
288 0 : sum=0.d0
289 0 : do j=max(1,i-mlmas),min(n,i+mumas)
290 0 : sum=sum+fmas(i-j+mbdiag,j)*ynew(j)
291 : end do
292 0 : ak(i)=ak(i)+sum
293 : end do
294 : end if
295 0 : call decsol(1,n,lde,e_1D,n,n,ak,ip,lrd,rpar_decsol,lid,ipar_decsol,ier)
296 0 : return
297 : !
298 : ! -----------------------------------------------------------
299 : !
300 : 13 continue
301 : ! --- b is a banded matrix, jacobian a full matrix, second order
302 0 : write(*,*) 'SLVROD: second order not supported'
303 0 : call mesa_error(__FILE__,__LINE__)
304 : ! if (not_stage1) then
305 : ! do i=1,m1
306 : ! ak(i)=ak(i)+ynew(i)
307 : ! end do
308 : ! do i=1,nm1
309 : ! sum=0.d0
310 : ! do j=max(1,i-mlmas),min(nm1,i+mumas)
311 : ! sum=sum+fmas(i-j+mbdiag,j)*ynew(j+m1)
312 : ! end do
313 : ! im1=i+m1
314 : ! ak(im1)=ak(im1)+sum
315 : ! end do
316 : ! end if
317 : ! if (ijob == 14) GOTO 45
318 : ! GOTO 48
319 : !
320 : ! -----------------------------------------------------------
321 : !
322 : 4 continue
323 : ! --- b is a banded matrix, jacobian a banded matrix
324 0 : if (not_stage1) then
325 0 : do i=1,n
326 0 : sum=0.d0
327 0 : do j=max(1,i-mlmas),min(n,i+mumas)
328 0 : sum=sum+fmas(i-j+mbdiag,j)*ynew(j)
329 : end do
330 0 : ak(i)=ak(i)+sum
331 : end do
332 : end if
333 0 : call decsol(1,n,lde,e_1D,mle,mue,ak,ip,lrd,rpar_decsol,lid,ipar_decsol,ier)
334 0 : return
335 : !
336 : ! -----------------------------------------------------------
337 : !
338 : 5 continue
339 : ! --- b is a full matrix, jacobian a full matrix
340 0 : if (not_stage1) then
341 0 : do i=1,n
342 : sum=0.d0
343 0 : do j=1,n
344 0 : sum=sum+fmas(i,j)*ynew(j)
345 : end do
346 0 : ak(i)=ak(i)+sum
347 : end do
348 : end if
349 0 : call decsol(1,n,lde,e_1D,n,n,ak,ip,lrd,rpar_decsol,lid,ipar_decsol,ier)
350 0 : return
351 : !
352 : ! -----------------------------------------------------------
353 : !
354 : 15 continue
355 : ! --- b is a full matrix, jacobian a full matrix, second order
356 0 : write(*,*) 'SLVROD: second order not supported'
357 0 : call mesa_error(__FILE__,__LINE__)
358 : ! if (not_stage1) then
359 : ! do i=1,m1
360 : ! ak(i)=ak(i)+ynew(i)
361 : ! end do
362 : ! do i=1,nm1
363 : ! sum=0.d0
364 : ! do j=1,nm1
365 : ! sum=sum+fmas(i,j)*ynew(j+m1)
366 : ! end do
367 : ! im1=i+m1
368 : ! ak(im1)=ak(im1)+sum
369 : ! end do
370 : ! end if
371 : ! GOTO 48
372 : !
373 : ! -----------------------------------------------------------
374 : !
375 : 6 continue
376 : ! --- b is a full matrix, jacobian a banded matrix
377 : ! --- this option is not provided
378 0 : if (not_stage1) then
379 0 : do i=1,n
380 : sum=0.d0
381 0 : do j=1,n
382 0 : sum=sum+fmas(i,j)*ynew(j)
383 : end do
384 0 : ak(i)=ak(i)+sum
385 : end do
386 0 : call decsol(1,n,lde,e_1D,mle,mue,ak,ip,lrd,rpar_decsol,lid,ipar_decsol,ier)
387 : end if
388 : return
389 : !
390 : ! -----------------------------------------------------------
391 : !
392 : 55 continue
393 0 : write(*,*) 'slvrod: invalid ijob', ijob
394 0 : call mesa_error(__FILE__,__LINE__) ! slvrod
395 81172 : end subroutine slvrod
|
