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