Line data Source code
1 0 : subroutine decomc(n,fjac,ldjac,fmas,ldmas,mlmas,mumas,
2 0 : & m1,m2,nm1,alphn,betan,e2r,e2i,lde1,ip2,br,bi,ierr,ijob,
3 : & mle,mue,mbjac,mbb,mdiag,mdiff,mbdiag,
4 0 : & decsolc,decsolcs,sparse_jac,nzmax,isparse,ia,ja,sar,sai,
5 : & lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol)
6 : use mtx_lib,only:find_loc_in_sparse
7 : implicit none
8 : interface
9 : #include "mtx_decsolc.dek"
10 : #include "mtx_decsolcs.dek"
11 : end interface
12 : integer :: m1, m2, nm1, lde1, ijob
13 : integer :: mle, mue, mbjac, mbb, mdiag, mdiff, mbdiag
14 : integer :: nzmax, isparse, lcd, lrd, lid
15 : integer :: ierr, ip2(nm1), n, ldjac, ldmas, mlmas, mumas
16 : integer :: ia(*), ja(nzmax) ! ia(n+1) when used; ia(2) when not.
17 : double precision :: sparse_jac(nzmax)
18 : double precision :: sar(nzmax), sai(nzmax)
19 : complex(dp), intent(inout), pointer :: cpar_decsol(:) ! (lcd)
20 : real(dp), intent(inout), pointer :: rpar_decsol(:) ! (lrd)
21 : integer, intent(inout), pointer :: ipar_decsol(:) ! (lid)
22 : double precision :: fjac(ldjac,n), fmas(ldmas,nm1)
23 : double precision :: e2r(lde1,nm1), e2i(lde1,nm1)
24 : double precision :: br(n), bi(n), alphn, betan
25 :
26 : ! LOCALS
27 : integer :: i, j, k, jm1, mm, imle, ib, hint
28 0 : double precision :: abno, alp, bet, sumr, sumi, sums, bb, ffma
29 :
30 0 : ierr = 0
31 0 : GOTO (1,2,3,4,5,6,55,8,9,55,11,12,13,14,15), ijob
32 : !
33 : ! -----------------------------------------------------------
34 : !
35 : 1 continue
36 : ! --- b=identity, jacobian a full matrix
37 0 : do j=1,n
38 0 : do i=1,n
39 0 : e2r(i,j)=-fjac(i,j)
40 0 : e2i(i,j)=0.d0
41 : end do
42 0 : e2r(j,j)=e2r(j,j)+alphn
43 0 : e2i(j,j)=betan
44 : end do
45 0 : call decsolc(0,n,lde1,e2r,e2i,n,n,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
46 0 : return
47 : !
48 : ! -----------------------------------------------------------
49 : !
50 : 8 continue
51 : ! --- b=identity, jacobian a sparse matrix
52 0 : do j=1,nzmax
53 0 : sar(j) = -sparse_jac(j)
54 : end do
55 0 : sai(1:nzmax) = 0
56 0 : do j=1,n
57 0 : do k=ia(j),ia(j+1)-1
58 0 : i = ja(k)
59 0 : if (i == j) then
60 0 : sar(k) = sar(k) + alphn
61 0 : sai(k) = betan
62 0 : exit
63 : end if
64 : end do
65 : end do
66 0 : call decsolcs(0,n,nzmax,ia,ja,sar,sai,br,bi,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
67 0 : return
68 : !
69 : ! -----------------------------------------------------------
70 : !
71 : 11 continue
72 : ! --- b=identity, jacobian a full matrix, second order
73 0 : do j=1,nm1
74 0 : jm1=j+m1
75 0 : do i=1,nm1
76 0 : e2r(i,j)=-fjac(i,jm1)
77 0 : e2i(i,j)=0.d0
78 : end do
79 0 : e2r(j,j)=e2r(j,j)+alphn
80 0 : e2i(j,j)=betan
81 : end do
82 0 : 45 mm=m1/m2
83 0 : abno=alphn**2+betan**2
84 0 : alp=alphn/abno
85 0 : bet=betan/abno
86 0 : do j=1,m2
87 0 : do i=1,nm1
88 0 : sumr=0.d0
89 0 : sumi=0.d0
90 0 : do k=0,mm-1
91 0 : sums=sumr+fjac(i,j+k*m2)
92 0 : sumr=sums*alp+sumi*bet
93 0 : sumi=sumi*alp-sums*bet
94 : end do
95 0 : e2r(i,j)=e2r(i,j)-sumr
96 0 : e2i(i,j)=e2i(i,j)-sumi
97 : end do
98 : end do
99 0 : call decsolc(0,nm1,lde1,e2r,e2i,nm1,nm1,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
100 0 : return
101 : !
102 : ! -----------------------------------------------------------
103 : !
104 : 2 continue
105 : ! --- b=identity, jacobian a banded matrix
106 0 : do j=1,n
107 0 : do i=1,mbjac
108 0 : imle=i+mle
109 0 : e2r(imle,j)=-fjac(i,j)
110 0 : e2i(imle,j)=0.d0
111 : end do
112 0 : e2r(mdiag,j)=e2r(mdiag,j)+alphn
113 0 : e2i(mdiag,j)=betan
114 : end do
115 0 : call decsolc(0,n,lde1,e2r,e2i,mle,mue,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
116 0 : return
117 : !
118 : ! -----------------------------------------------------------
119 : !
120 : 12 continue
121 : ! --- b=identity, jacobian a banded matrix, second order
122 0 : do j=1,nm1
123 0 : jm1=j+m1
124 0 : do i=1,mbjac
125 0 : e2r(i+mle,j)=-fjac(i,jm1)
126 0 : e2i(i+mle,j)=0.d0
127 : end do
128 0 : e2r(mdiag,j)=e2r(mdiag,j)+alphn
129 0 : e2i(mdiag,j)=e2i(mdiag,j)+betan
130 : end do
131 0 : 46 mm=m1/m2
132 0 : abno=alphn**2+betan**2
133 0 : alp=alphn/abno
134 0 : bet=betan/abno
135 0 : do j=1,m2
136 0 : do i=1,mbjac
137 0 : sumr=0.d0
138 0 : sumi=0.d0
139 0 : do k=0,mm-1
140 0 : sums=sumr+fjac(i,j+k*m2)
141 0 : sumr=sums*alp+sumi*bet
142 0 : sumi=sumi*alp-sums*bet
143 : end do
144 0 : imle=i+mle
145 0 : e2r(imle,j)=e2r(imle,j)-sumr
146 0 : e2i(imle,j)=e2i(imle,j)-sumi
147 : end do
148 : end do
149 0 : call decsolc(0,nm1,lde1,e2r,e2i,mle,mue,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
150 0 : return
151 : !
152 : ! -----------------------------------------------------------
153 : !
154 : 3 continue
155 : ! --- b is a banded matrix, jacobian a full matrix
156 0 : do j=1,n
157 0 : do i=1,n
158 0 : e2r(i,j)=-fjac(i,j)
159 0 : e2i(i,j)=0.d0
160 : end do
161 : end do
162 0 : do j=1,n
163 0 : do i=max(1,j-mumas),min(n,j+mlmas)
164 0 : bb=fmas(i-j+mbdiag,j)
165 0 : e2r(i,j)=e2r(i,j)+alphn*bb
166 0 : e2i(i,j)=betan*bb
167 : end do
168 : end do
169 0 : call decsolc(0,n,lde1,e2r,e2i,n,n,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
170 0 : return
171 : !
172 : ! -----------------------------------------------------------
173 : !
174 : 9 continue
175 : ! --- b is a banded matrix, jacobian a sparse matrix
176 0 : do j=1,nzmax
177 0 : sar(j) = -sparse_jac(j)
178 : end do
179 0 : sai(1:nzmax) = 0
180 0 : do j=1,n
181 0 : hint = 0
182 0 : do i=max(1,j-mumas),min(n,j+mlmas)
183 : ! set k = location in sa for (i,j)
184 : ierr = 0
185 0 : call find_loc_in_sparse(isparse,n,nzmax,ia,ja,i,j,hint,k,ierr)
186 0 : if (ierr == 0) then
187 0 : bb=fmas(i-j+mbdiag,j)
188 0 : sar(k)=sar(k)+alphn*bb
189 0 : sai(k)=betan*bb
190 0 : else if (fmas(i-j+mbdiag,j) /= 0) then
191 : return ! bad sparsity
192 : else
193 : ierr = 0
194 : end if
195 0 : hint = k
196 : end do
197 : end do
198 : !write(*,*) 'decomc call decsolcs'
199 0 : call decsolcs(0,n,nzmax,ia,ja,sar,sai,br,bi,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
200 : !write(*,*) 'decomc back decsolcs'
201 0 : return
202 : !
203 : ! -----------------------------------------------------------
204 : !
205 : 13 continue
206 : ! --- b is a banded matrix, jacobian a full matrix, second order
207 0 : do j=1,nm1
208 0 : jm1=j+m1
209 0 : do i=1,nm1
210 0 : e2r(i,j)=-fjac(i,jm1)
211 0 : e2i(i,j)=0.d0
212 : end do
213 0 : do i=max(1,j-mumas),min(nm1,j+mlmas)
214 0 : ffma=fmas(i-j+mbdiag,j)
215 0 : e2r(i,j)=e2r(i,j)+alphn*ffma
216 0 : e2i(i,j)=e2i(i,j)+betan*ffma
217 : end do
218 : end do
219 : GOTO 45
220 : !
221 : ! -----------------------------------------------------------
222 : !
223 : 4 continue
224 : ! --- b is a banded matrix, jacobian a banded matrix
225 0 : do j=1,n
226 0 : do i=1,mbjac
227 0 : imle=i+mle
228 0 : e2r(imle,j)=-fjac(i,j)
229 0 : e2i(imle,j)=0.d0
230 : end do
231 0 : do i=max(1,mumas+2-j),min(mbb,mumas+1-j+n)
232 0 : ib=i+mdiff
233 0 : bb=fmas(i,j)
234 0 : e2r(ib,j)=e2r(ib,j)+alphn*bb
235 0 : e2i(ib,j)=betan*bb
236 : end do
237 : end do
238 0 : call decsolc(0,n,lde1,e2r,e2i,mle,mue,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
239 0 : return
240 : !
241 : ! -----------------------------------------------------------
242 : !
243 : 14 continue
244 : ! --- b is a banded matrix, jacobian a banded matrix, second order
245 0 : do j=1,nm1
246 0 : jm1=j+m1
247 0 : do i=1,mbjac
248 0 : e2r(i+mle,j)=-fjac(i,jm1)
249 0 : e2i(i+mle,j)=0.d0
250 : end do
251 0 : do i=1,mbb
252 0 : ib=i+mdiff
253 0 : ffma=fmas(i,j)
254 0 : e2r(ib,j)=e2r(ib,j)+alphn*ffma
255 0 : e2i(ib,j)=e2i(ib,j)+betan*ffma
256 : end do
257 : end do
258 : GOTO 46
259 : !
260 : ! -----------------------------------------------------------
261 : !
262 : 5 continue
263 : ! --- b is a full matrix, jacobian a full matrix
264 0 : do j=1,n
265 0 : do i=1,n
266 0 : bb=fmas(i,j)
267 0 : e2r(i,j)=bb*alphn-fjac(i,j)
268 0 : e2i(i,j)=bb*betan
269 : end do
270 : end do
271 0 : call decsolc(0,n,lde1,e2r,e2i,n,n,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
272 0 : return
273 : !
274 : ! -----------------------------------------------------------
275 : !
276 : 15 continue
277 : ! --- b is a full matrix, jacobian a full matrix, second order
278 0 : do j=1,nm1
279 0 : jm1=j+m1
280 0 : do i=1,nm1
281 0 : e2r(i,j)=alphn*fmas(i,j)-fjac(i,jm1)
282 0 : e2i(i,j)=betan*fmas(i,j)
283 : end do
284 : end do
285 : GOTO 45
286 : !
287 : ! -----------------------------------------------------------
288 : !
289 : 6 continue
290 : ! --- b is a full matrix, jacobian a banded matrix
291 : ! --- this option is not provided
292 : return
293 : !
294 : ! -----------------------------------------------------------
295 : !
296 : 55 continue
297 0 : write(*,*) 'decomc: invalid ijob', ijob
298 0 : call mesa_error(__FILE__,__LINE__) ! decomc
299 : end subroutine decomc
300 : !
301 : ! ***********************************************************
302 :
303 :
304 0 : subroutine decsolc_done(n,fjac,ldjac,fmas,ldmas,mlmas,mumas,
305 0 : & m1,m2,nm1,alphn,betan,e2r,e2i,lde1,ip2,br,bi,ierr,ijob,
306 : & mle,mue,mbjac,mbb,mdiag,mdiff,mbdiag,
307 0 : & decsolc,decsolcs,sparse_jac,nzmax,isparse,ia,ja,sar,sai,
308 : & lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol)
309 : implicit none
310 : interface
311 : #include "mtx_decsolc.dek"
312 : #include "mtx_decsolcs.dek"
313 : end interface
314 : integer :: m1, m2, nm1, lde1, ijob
315 : integer :: mle, mue, mbjac, mbb, mdiag, mdiff, mbdiag
316 : integer :: nzmax, isparse, lcd, lrd, lid
317 : integer :: ierr, ip2(nm1), n, ldjac, ldmas, mlmas, mumas
318 : integer :: ia(:) ! (n+1)
319 : integer :: ja(:) ! (nzmax)
320 : real(dp) :: sparse_jac(:) ! (nzmax)
321 : real(dp) :: sar(:) ! (nzmax)
322 : real(dp) :: sai(:) ! (nzmax)
323 : complex(dp), intent(inout), pointer :: cpar_decsol(:) ! (lcd)
324 : real(dp), intent(inout), pointer :: rpar_decsol(:) ! (lrd)
325 : integer, intent(inout), pointer :: ipar_decsol(:) ! (lid)
326 :
327 : double precision :: fjac(ldjac,n), fmas(ldmas,nm1)
328 : double precision :: e2r(lde1,nm1), e2i(lde1,nm1)
329 : double precision :: br(n), bi(n), alphn, betan
330 :
331 :
332 0 : GOTO (1,2,3,4,5,6,55,8,9,55,11,12,13,14,15), ijob
333 : !
334 : ! -----------------------------------------------------------
335 : !
336 : 1 continue
337 : ! --- b=identity, jacobian a full matrix
338 0 : call decsolc(2,n,lde1,e2r,e2i,n,n,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
339 0 : return
340 : !
341 : ! -----------------------------------------------------------
342 : !
343 : 8 continue
344 : ! --- b=identity, jacobian a sparse matrix
345 0 : call decsolcs(2,n,nzmax,ia,ja,sar,sai,br,bi,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
346 0 : return
347 : !
348 : ! -----------------------------------------------------------
349 : !
350 : 11 continue
351 : ! --- b=identity, jacobian a full matrix, second order
352 : 45 continue
353 0 : call decsolc(2,nm1,lde1,e2r,e2i,nm1,nm1,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
354 0 : return
355 : !
356 : ! -----------------------------------------------------------
357 : !
358 : 2 continue
359 : ! --- b=identity, jacobian a banded matrix
360 0 : call decsolc(2,n,lde1,e2r,e2i,mle,mue,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
361 0 : return
362 : !
363 : ! -----------------------------------------------------------
364 : !
365 : 12 continue
366 : ! --- b=identity, jacobian a banded matrix, second order
367 : 46 continue
368 0 : call decsolc(2,nm1,lde1,e2r,e2i,mle,mue,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
369 0 : return
370 : !
371 : ! -----------------------------------------------------------
372 : !
373 : 3 continue
374 : ! --- b is a banded matrix, jacobian a full matrix
375 0 : call decsolc(2,n,lde1,e2r,e2i,n,n,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
376 0 : return
377 : !
378 : ! -----------------------------------------------------------
379 : !
380 : 9 continue
381 : ! --- b is a banded matrix, jacobian a sparse matrix
382 0 : call decsolcs(2,n,nzmax,ia,ja,sar,sai,br,bi,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
383 0 : return
384 : !
385 : ! -----------------------------------------------------------
386 : !
387 : 13 continue
388 : ! --- b is a banded matrix, jacobian a full matrix, second order
389 : GOTO 45
390 : !
391 : ! -----------------------------------------------------------
392 : !
393 : 4 continue
394 : ! --- b is a banded matrix, jacobian a banded matrix
395 0 : call decsolc(2,n,lde1,e2r,e2i,mle,mue,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
396 0 : return
397 : !
398 : ! -----------------------------------------------------------
399 : !
400 : 14 continue
401 : ! --- b is a banded matrix, jacobian a banded matrix, second order
402 : GOTO 46
403 : !
404 : ! -----------------------------------------------------------
405 : !
406 : 5 continue
407 : ! --- b is a full matrix, jacobian a full matrix
408 0 : call decsolc(2,n,lde1,e2r,e2i,n,n,br,bi,ip2,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ierr)
409 0 : return
410 : !
411 : ! -----------------------------------------------------------
412 : !
413 : 15 continue
414 : ! --- b is a full matrix, jacobian a full matrix, second order
415 : GOTO 45
416 : !
417 : ! -----------------------------------------------------------
418 : !
419 : 6 continue
420 : ! --- b is a full matrix, jacobian a banded matrix
421 : ! --- this option is not provided
422 : return
423 : !
424 : ! -----------------------------------------------------------
425 : !
426 : 55 continue
427 0 : write(*,*) 'decsolc_done: invalid ijob', ijob
428 0 : call mesa_error(__FILE__,__LINE__) ! decsolc_done
429 : end subroutine decsolc_done
|
