ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
spider_solver.f
Go to the documentation of this file.
1  !subroutine solve_fft(isol,wdm)
2  subroutine solve(isol,wdm)
3  include 'double.inc'
4  include 'param.inc'
5  include 'comblc.inc'
6  common /comffs/ axm(nip),ax0(nip),axp(nip)
7 
8  real time_beg,time_end
9  real*8 wdm(nip,njp)
10  real*8 wdmx(nip,njp)
11  dimension fwrk(nip,njp),ffwrk(nip,njp)
12 
13  do i=1,ni
14  do j=1,nj
15  wdmx(i,j)=wdm(i,j)
16  enddo
17  enddo
18 
19  do i=1,ni
20  do j=1,nj
21  fwrk(i,j)=-curf(i,j)*r(i)
22  enddo
23  enddo
24 
25  call foursol(ni,nj,wdm,fwrk,ffwrk,r,z,axm,ax0,axp)
26 
27  wght=0.5d0
28  if(itin.le.2 ) wght=1.0d0
29 
30  do i=2,ni1
31  do j=2,nj1
32  wdm(i,j)=wdm(i,j)*wght+(1.d0-wght)*wdmx(i,j)
33  enddo
34  enddo
35 
36  return
37  end
38 
39 !*******************************************************************
40 
41  !subroutine solve(isol,wdm)
42  subroutine solve_spar(isol,wdm)
43 
44  include 'double.inc'
45  include 'param.inc'
46  include 'comblc.inc'
47 
48  common /comwrc/ rsp,p,ip
49 
50  real*8 zw(neqp),rsp(nspp),wdm(nip,njp)
51 
52  integer p(neqp),ip(neqp),isp(nspp),ipath,flag,esp
53 
54  equivalence(rsp(1),isp(1))
55 
56 c equivalence(right(1),zw(1))
57 
58  ipath=3
59 
60  if(isol.ne.0) go to 20
61 
62  call odrvd(neq,ia,ja,a,p,ip,nspp,isp,1,flag)
63 
64 c do 10 i=1,neqp
65 c ip(i)=i
66 c p(i)=i
67 c10 continue
68 
69  !write(6,*) 'odrv flag=',flag
70 
71  ipath=1
72  isol=1
73 
74  20 continue
75 
76 
77  call sdrvd(neq,p,ip,ia,ja,a,right,zw,nspp,
78  * isp,rsp,esp,ipath,flag)
79 
80 c do 860 i=1,neq
81 c write(6,*) 'zw(i) i',i,zw(i)
82 c860 continue
83 c
84  !write(6,*) 'sdrv: flag,esp',flag,esp
85 c
86 cccc call nev(zw)
87 
88 c raspakovka reshenia
89 
90  wght=0.5d0
91  if(itin.le.2 ) wght=1.0d0
92 
93  do 500 i=2,ni1
94  do 500 j=2,nj1
95 
96  ieq=nlin(i,j)
97 
98  wdm(i,j)=zw(ieq)*wght+(1.d0-wght)*wdm(i,j)
99 
100  500 continue
101 
102  return
103  end
104 
105  subroutine nev(zw)
106 
107  include 'double.inc'
108  include 'param.inc'
109  include 'comblc.inc'
110 
111  real*8 zw(neqp)
112 
113  do 10 il=1,neq
114 
115  i1=ia(il)
116  i2=ia(il+1)-1
117 
118  znev=0.d0
119  znmx=0.d0
120 
121  do 100 im=i1,i2
122 
123  ic=ja(im)
124 
125  znev=znev+a(im)*zw(ic)
126 
127  ! if(il.eq.1) then
128  !
129  ! write(6,*) 'im a(im)',im,a(im)
130  ! write(6,*) 'ic zw(ic)',ic,zw(ic)
131  ! write(6,*) 'znev ' ,znev
132  ! endif
133 
134  100 continue
135 
136  znev=right(il)-znev
137  znab=dabs(znev)
138 c write(6,*) '***right znev ***',il,right(il),znev
139  znmx=dmax1(znmx,znab)
140 
141  10 continue
142 
143 c write(6,*) 'nev:',znmx
144 
145  return
146  end
147 c-----------------------------------------------------------------------
148 !***************************************************************
149 
150  subroutine foursol(ni,nj,u,f,ff,x,y,axm,ax0,axp)
151 
152  implicit real*8(a-h,o-z)
153  include 'param.inc'
154 
155  dimension u(nip,njp),f(nip,njp),x(nip),y(njp)
156  dimension ff(nip,njp)
157  dimension ff1(njp),w(njp),slam(njp),a(njp)
158  dimension ap(nip),bp(nip),cp(nip),fp(nip)
159  dimension alf(nip),bet(nip),u1(nip)
160  dimension axm(nip),axp(nip),ax0(nip)
161 
162 c-------power nj=2*mm definition
163  md=nj-1
164  do m=1,100
165  mm=m
166  if(md/2*2.ne.md) then
167  write(*,*) 'nj-1 must be power of 2.program is terminated '
168  write(*,*) 'nj=', nj
169  stop
170  endif
171  md=md/2
172  if(md.eq.1) go to 446
173  enddo
174  446 continue
175  nj1=nj-1
176  ni1=ni-1
177 
178  hy=y(2)-y(1)
179 
180  do i=2,ni1
181  do j=2,nj1
182  u(i,j)=f(i,j)
183  enddo
184  enddo
185 
186  do i=2,ni-1
187  u(i,2)=u(i,2)-u(i,1)/hy**2
188  u(i,nj-1)=u(i,nj-1)-u(i,nj)/hy**2
189  enddo
190 
191  do j=2,nj-1
192  u(2,j)=u(2,j)-u(1,j)*axm(2)
193  u(ni-1,j)=u(ni-1,j)-u(ni,j)*axp(ni1)
194  enddo
195 !*************************************************
196  do i=2,ni1
197 
198  do j=2,nj1
199  w(j)=u(i,j)
200  enddo
201 
202  call rft1(mm,w,a,3)
203 
204  do k=2,nj1
205  ff(i,k)=a(k)
206  enddo
207 
208  enddo
209 
210 !*************************************************
211 
212  do k=2,nj-1
213 
214  slam(k)=4.d0*dsin(0.5d0*(k-1.d0)*pi/nj1)**2/hy**2
215 
216  do i=2,ni-1
217  ap(i-1)=axp(i)
218  bp(i-1)=axm(i)
219  cp(i-1)=ax0(i)+slam(k)
220  fp(i-1)=-ff(i,k)
221  enddo
222 
223  ap(ni-2)=0.d0
224  bp(1)=0.d0
225 
226  call prog1d(ni-2,u1(2), ap,bp,cp,fp, alf,bet)
227 
228  do i=2,ni-1
229  ff(i,k)=u1(i)
230  enddo
231 
232  enddo
233 !*************************************************
234  do i=2,ni-1
235  do k=2,nj-1
236  a(k)=ff(i,k)
237  enddo
238 
239  call rft1(mm,a,w,6)
240 
241  do j=2,nj-1
242  u(i,j)=w(j)
243  enddo
244  enddo
245 
246 ! zerr=0.d0
247 ! do i=2,ni-1
248 ! do j=2,nj-1
249 ! znew=axm(i)*u(i-1,j)-ax0(i)*u(i,j)+axp(i)*u(i+1,j)
250 ! + +(u(i,j-1)-2.d0*u(i,j)+u(i,j+1))/hy**2-f(i,j)
251 ! ernew=dabs(znew)
252 ! ernew=dmax1(zerr,ernew)
253 ! enddo
254 ! enddo
255 !
256 ! write(6,*) 'fursol:new',ernew
257 
258  return
259  end
260 !**********************************************************************
261 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
262 
263  SUBROUTINE rft1 (MM, X, Y, MODE)
264 c fast fourier transform in one dimension.
265 c ==================================================================
266  IMPLICIT REAL*8(a-h,o-z)
267  dimension x(*), y(*)
268  COMMON /fdata/f,rttwo, n, n2, m, nc, ns, nd, nr
269  COMMON /fwork/ w(2048)
270 c ==================================================================
271 c==== fixed DATA preparation.
272  md = mm
273  IF(mode .NE. 1 .AND. mode .NE. 4) md = mm + 1
274  CALL setft(md)
275  n1 = 2 ** mm + 1
276 c==== copy input vector into working storage.
277  DO 10 i = 1, n1
278  10 w(nd+i) = x(i) * f
279  go to(20, 30, 40, 50, 30, 40), mode
280 c==== general analysis.
281  20 CALL rpa
282  w(nr+n1) = 0.0d0
283  go to 60
284 c==== dcodsine analysis or synthesis.
285  30 CALL rca
286  go to 60
287 c==== dsine analysis or synthesis.
288  40 CALL rsa
289  w(nr+1) = x(1)
290  w(nr+n1) = x(n1)
291  go to 60
292 c==== general synthesis.
293  50 CALL rps
294  w(nr+n1) = w(nr+1)
295 c==== copy into output vector.
296  60 DO 70 i = 1, n1
297  70 y(i) = w(nr+i)
298  RETURN
299 c ==================================================================
300  END
301  SUBROUTINE rpa
302 c REAL analysis.
303 c ==================================================================
304  IMPLICIT REAL*8(a-h,o-z)
305  COMMON /fdata/ f,rttwo,n, n2, m, nc, ns, ndata, nres
306  COMMON /fwork/ w(2048)
307 c ==================================================================
308  nq = n2
309  nir = ndata
310  nor = ndata + n
311  DO 100 l = 1, m
312  nii = nir + n2
313  noi = nor + n2
314 c==== tprime = 0.
315  DO 10 it = 1, nq
316  w(nor+it) = w(nir+it) + w(nir+it+nq)
317  10 w(noi+it) = w(nir+it) - w(nir+it+nq)
318 c==== 0 .LT. tprime .LE. np / 2
319  no1 = nq
320  no2 = n2 - nq
321  IF(no1 - no2) 20, 40, 60
322  20 ni1 = no1 + no1
323  ni2 = ni1 + nq
324  cc = w(nc+no1)
325  ss = w(ns+no1)
326  DO 30 it = 1, nq
327  re = cc * w(nir+ni2+it) - ss * w(nii+ni2+it)
328  ai = ss * w(nir+ni2+it) + cc * w(nii+ni2+it)
329  w(nor+no1+it) = w(nir+ni1+it) + re
330  w(nor+no2+it) = w(nir+ni1+it) - re
331  w(noi+no1+it) = w(nii+ni1+it) + ai
332  30 w(noi+no2+it) = - w(nii+ni1+it) + ai
333  no1 = no1 + nq
334  no2 = no2 - nq
335  IF(no1 - no2) 20, 40, 60
336 c==== tprime = np/2.
337  40 DO 50 it = 1, nq
338  w(nor+no1+it) = w(nii +it)
339  50 w(noi+no1+it) = w(nii+nq+it)
340  60 nt = nir
341  nir = nor
342  nor = nt
343  nq = nq / 2
344  100 CONTINUE
345  nres = nir
346  RETURN
347 c ==================================================================
348  END
349  SUBROUTINE rps
350 c REAL synthesis.
351 c ==================================================================
352  IMPLICIT REAL*8(a-h,o-z)
353  COMMON /fdata/ f,rttwo,n, n2, m, nc, ns, ndata, nres
354  COMMON /fwork/ w(2048)
355 c ==================================================================
356  nq = 1
357  nir = ndata
358  nor = ndata + n
359  l = m
360  10 nii = nir + n2
361  noi = nor + n2
362 c==== tprime = 0.
363  DO 20 it = 1, nq
364  w(nor+it) = w(nir+it) + w(nii+it)
365  20 w(nor+it+nq) = w(nir+it) - w(nii+it)
366 c==== 0 .LT. tprime .LT. np/2.
367  ni1 = nq
368  ni2 = n2 - nq
369  IF(ni1 - ni2) 40, 80, 120
370  40 no1 = ni1 + ni1
371  no2 = no1 + nq
372  cc = w(nc+ni1)
373  ss = w(ns+ni1)
374  DO 50 it = 1, nq
375  w(nor+no1+it) = w(nir+ni1+it) + w(nir+ni2+it)
376  re = w(nir+ni1+it) - w(nir+ni2+it)
377  w(noi+no1+it) = w(nii+ni1+it) - w(nii+ni2+it)
378  ai = w(nii+ni1+it) + w(nii+ni2+it)
379  w(nor+no2+it) = cc * re + ss * ai
380  50 w(noi+no2+it) = cc * ai - ss * re
381  ni1 = ni1 + nq
382  ni2 = ni2 - nq
383  IF(ni1 - ni2) 40, 80, 120
384  80 DO 90 it = 1, nq
385  w(noi +it) = w(nir+ni1+it) * 2.0d0
386  90 w(noi+nq+it) = w(nii+ni1+it) * 2.0d0
387  120 nt = nir
388  nir = nor
389  nor = nt
390  nq = nq + nq
391  l = l - 1
392  IF(l .GT. 0) go to 10
393  nres = nir
394  RETURN
395 c ==================================================================
396  END
397  SUBROUTINE rsa
398 c REAL odd analysis or synthesis.
399 c ==================================================================
400  IMPLICIT REAL*8(a-h,o-z)
401  COMMON /fdata/ f,rttwo,n, n2, m, nc, ns, ndata, nres
402  COMMON /fwork/ w(2048)
403 c ==================================================================
404  nq = n2
405  nib = ndata + 1
406  nob = ndata + n2 + 2
407  DO 200 l = 1, m
408  nim = nib + nq
409  nie = nim + nq
410  nqh = nq / 2
411  nql = (nq - 1) / 2
412  nom = nob + nql + 1
413 c==== tprime = 0.
414  IF(nql .EQ. 0) go to 20
415  DO 10 it = 1, nql
416  w(nob+it) = w(nib+it) - w(nim-it)
417  10 w(nom+it) = w(nib+it) + w(nim-it)
418  20 IF(nqh .NE. nql) w(nom) = 2.0d0 * w(nib+nqh)
419 c==== 0 .LT. tprime .LT. np/2.
420  no1 = nq
421  no2 = n2 - nq
422  IF(no1 - no2) 40, 80, 120
423  40 ni = no1 + no1
424  w(nob+no1) = w(nib+ni) + w(nim+ni)
425  w(nob+no2) = - w(nib+ni) + w(nim+ni)
426  IF(nql .EQ. 0) go to 60
427  cc = w(nc+no1)
428  ss = w(ns+no1)
429  DO 50 it=1,nql
430  re = cc * w(nie+ni-it) - ss * w(nim+ni-it)
431  ai = ss * w(nie+ni-it) + cc * w(nim+ni-it)
432  w(nom+no1+it) = w(nim+ni+it) - re
433  w(nom+no2+it) = w(nim+ni+it) + re
434  w(nob+no1+it) = w(nib+ni+it) + ai
435  50 w(nob+no2+it) = - w(nib+ni+it) + ai
436  60 IF(nqh .EQ. nql) go to 70
437  nr = no1 / 2
438  w(nom+no1)=2.d0*(w(ns+nr)*w(nim+ni+nqh) + w(nc+nr)*w(nib+ni+nqh))
439  w(nom+no2)=2.d0*(w(nc+nr)*w(nim+ni+nqh) - w(ns+nr)*w(nib+ni+nqh))
440  70 no1 = no1 + nq
441  no2 = no2 - nq
442  IF(no1 - no2) 40, 80, 120
443  80 w(nob+no1) = w(nim)
444  IF(nql .EQ. 0) go to 100
445  DO 90 it = 1, nql
446  w(nom+no1+it) = w(nim+it)
447  90 w(nob+no1+it) = w(nie-it)
448  100 IF(nqh .NE. nql) w(nom+no1) = rttwo * w(nim+nqh)
449  120 nt = nib
450  nib = nob
451  nob = nt
452  nq = nqh
453  200 CONTINUE
454  nres = nib - 1
455  RETURN
456 c ==================================================================
457  END
458  SUBROUTINE rca
459 c REAL even analysis or synthesis.
460 c ==================================================================
461  IMPLICIT REAL*8(a-h,o-z)
462  COMMON /fdata/ f,rttwo,n, n2, m, nc, ns, ndata, nres
463  COMMON /fwork/ w(2048)
464 c ==================================================================
465  nq = n2
466  nib = ndata + 1
467  nob = ndata + n2 + 2
468  DO 200 l = 1, m
469  nim = nib + nq
470  nie = nim + nq
471  nqh = nq / 2
472  nql = (nq - 1) / 2
473  nom = nob + nql + 1
474 c==== tprime = 0.
475  w(nob ) = w(nib) + w(nim)
476  w(nob+n2) = w(nib) - w(nim)
477  IF(nql .EQ. 0) go to 20
478  DO 10 it = 1, nql
479  w(nob+it) = w(nib+it) + w(nim-it)
480  10 w(nom+it) = w(nib+it) - w(nim-it)
481  20 IF(nqh .NE. nql) w(nom) = w(nib+nqh) * 2.d0
482 c==== 0 .LT. tprime .LT. np / 2
483  no1 = nq
484  no2 = n2 - nq
485  IF(no1 - no2) 40, 80, 120
486  40 ni = no1 + no1
487  w(nob+no1) = w(nib+ni) + w(nim+ni)
488  w(nob+no2) = w(nib+ni) - w(nim+ni)
489  IF(nql .EQ. 0) go to 60
490  cc = w(nc+no1)
491  ss = w(ns+no1)
492  DO 50 it = 1, nql
493  re = cc * w(nim+ni-it) - ss * w(nie+ni-it)
494  ai = ss * w(nim+ni-it) + cc * w(nie+ni-it)
495  w(nob+no1+it) = w(nib+ni+it) + re
496  w(nob+no2+it) = w(nib+ni+it) - re
497  w(nom+no1+it) = w(nim+ni+it) - ai
498  50 w(nom+no2+it) = - w(nim+ni+it) - ai
499  60 IF(nqh .EQ. nql) go to 70
500  nr = no1 / 2
501  w(nom+no1)=2.d0*(w(nc+nr)*w(nib+ni+nqh) - w(ns+nr)*w(nim+ni+nqh))
502  w(nom+no2)=2.d0*(w(ns+nr)*w(nib+ni+nqh) + w(nc+nr)*w(nim+ni+nqh))
503  70 no1 = no1 + nq
504  no2 = no2 - nq
505  IF(no1 - no2) 40, 80, 120
506 c==== tprime = np/2
507  80 w(nob+no1) = w(nib+n2)
508  IF(nql .EQ. 0) go to 100
509  DO 90 it = 1, nql
510  w(nob+no1+it) = w(nim+it)
511  90 w(nom+no1+it) = - w(nie-it)
512  100 IF(nqh .NE. nql) w(nom+no1) = rttwo * w(nim+nqh)
513  120 nt = nib
514  nib = nob
515  nob = nt
516  nq = nqh
517  200 CONTINUE
518  nres = nib - 1
519  RETURN
520 c ==================================================================
521  END
522  SUBROUTINE setft (MM)
523 c prepare fixed DATA for fast fourier transform.
524 c ==================================================================
525  IMPLICIT REAL*8(a-h,o-z)
526  COMMON /fdata/f,rttwo, n, n2, m, nc, ns, ndata, nres
527  COMMON /fwork/ w(2048)
528 c DATA (m = - 1)
529  m = - 1
530 c ==================================================================
531  IF(mm .LT. 0) go to 90
532  IF(mm .EQ. m) RETURN
533  m = mm
534  n = 2 ** m
535  n2 = n / 2
536  f = 1.d0 / dfloat(n)
537  da = 6.2831853071796d0 * f
538  rttwo = 2.0d0
539  rttwo = dsqrt(rttwo)
540  f = dsqrt(f)
541  n1 = (n2 - 1) / 2
542  nc = 0
543  ns = n1
544  ndata = n1 + n1
545  IF(n1 .LE. 0) RETURN
546  DO 10 i = 1, n1
547  a = dfloat(i) * da
548  w(i) = dcos(a)
549  10 w(n1+i) = dsin(a)
550  RETURN
551  90 CONTINUE
552 c ==================================================================
553  END
554 
555 
556 
557 
558 
nlin
function nlin(i, j)
Definition: Matrix.f:173
rft1
subroutine rft1(MM, X, Y, MODE)
Definition: spider_solver.f:263
foursol
subroutine foursol(ni, nj, u, f, ff, x, y, axm, ax0, axp)
Definition: spider_solver.f:150
solve
subroutine solve(isol, wdm)
Definition: spider_solver.f:2
solve_spar
subroutine solve_spar(isol, wdm)
Definition: spider_solver.f:42
md
subroutine md
Definition: SPAR.f:347
odrvd
subroutine odrvd
Definition: SPAR.f:9
prog1d
subroutine prog1d(M, U, A, B, C, F, ALF, BET)
Definition: scu.f:822
output
subroutine output(NGRID, betpol, zli3)
Definition: Eq2_m.f:1
rpa
subroutine rpa
Definition: spider_solver.f:301
nev
subroutine nev(zw)
Definition: spider_solver.f:105
setft
subroutine setft(MM)
Definition: spider_solver.f:522
rca
subroutine rca
Definition: spider_solver.f:458
rps
subroutine rps
Definition: spider_solver.f:349
phys_functions::p
real(r8) function p(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:383
rsa
subroutine rsa
Definition: spider_solver.f:397
sdrvd
subroutine sdrvd
Definition: SPAR.f:1411