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