ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/JSOC/proj/globalhs/apps/inv2d.f
Revision: 1.3
Committed: Tue Jan 17 18:43:40 2017 UTC (6 years, 8 months ago) by tplarson
Branch: MAIN
CVS Tags: globalhs_version_23, globalhs_version_22, globalhs_version_21, globalhs_version_20, Ver_LATEST, globalhs_version_24, Ver_9-41, Ver_9-5, globalhs_version_19, Ver_9-1, Ver_9-2, globalhs_version_16, globalhs_version_17, globalhs_version_15, globalhs_version_18, Ver_9-4, Ver_9-3, Ver_9-0, HEAD
Changes since 1.2: +5 -0 lines
Log Message:
update compilation and version tracking for inv2d

File Contents

# Content
1 program inv2d
2 c 070701 Changed output format to accommodate l=1000
3
4 implicit double precision (a-h,o-z)
5 include 'include1.f'
6 include 'include2.f'
7
8 parameter (maxpoints=maxr1*maxt1,nblock1=10)
9
10 double precision split(maxnlm),spliti(maxl)
11 double precision sigma(maxnlm),sigmai(maxl)
12 double precision sigma1(maxnlm),ss1(maxnlm),ss2(maxnlm)
13 double precision rmesh1(maxr1),tmesh1(maxt1)
14 double precision rweights(maxr1),tweights(maxt1)
15 double precision fr(maxr1),ft(maxr1),fwr(maxr1),fwt(maxt1)
16 character*80 cefunc,outfile,splitfile
17 real dummy(2),dtime,etime,time1
18 real mur,mut
19 real f1(maxmodes,maxr1),f2(maxmodes,maxr1)
20 real g1(maxlm,maxt1),g2(maxlm,maxt1)
21 double precision f1x(maxmodes),f2x(maxmodes)
22 double precision g1x(maxlm),g2x(maxlm)
23 double precision f1i(maxr1),f2i(maxr1),g1i(maxt1),g2i(maxt1)
24 double precision f1o(maxt1),f2o(maxt1)
25 double precision a(0:maxs1-1,maxpoints)
26 double precision ax(0:maxs1-1),bx(nblock,nblock)
27 integer iax(maxpoints)
28 double precision b(maxpoints,maxpoints),right(maxpoints)
29 double precision c1(maxr1,maxr1),c2(maxt1,maxt1)
30 double precision d(maxpoints,maxpoints),help(3*maxpoints)
31 integer iwork(maxpoints)
32 c double precision d1(maxpoints,maxpoints),cov(maxpoints,maxpoints)
33 double precision var(maxpoints),dib(maxpoints)
34 double precision dibx(maxpoints,nblock1)
35 c double precision di(maxpoints)
36 double precision coeff(maxnlm)
37 double precision omega(maxr1,maxt1)
38 double precision fro(maxr1,maxt1),fto(maxr1,maxt1)
39 integer jav(maxpoints),kav(maxpoints)
40 real avkern(maxtav,maxrav)
41 integer nnlm1(maxnlm),llist(maxnlm)
42 real fnlm(maxnlm)
43 logical found
44 integer ijk(maxr1,maxt1)
45 logical ibig(maxpoints)
46 c integer lsplit(maxnlm)
47 c integer msplit(maxnlm),isplit(maxnlm),lmsplit(maxnlm)
48 integer mindex(maxnlm)
49 real x1(maxtav),x2(maxrav),yy2(maxtav,maxrav),mjta(4)
50 real xr(10)
51 integer ixl(maxmodes),ixh(maxmodes)
52
53
54 c nmode,lmode: identification of modes
55 c freq: frequency of mode (muHz)
56 c radmesh: mesh in radius for kernels and eigenfunctions
57 c nexp: number of expansion coefficients
58 c cefunc: control file for reading of eigenfunctions
59 c iomega: 0 to read omega from input
60 c 1 to read omega from file
61 c omegafile: file from which omega is read
62 c nrin, ntin: number of points where omega is given
63 c rin,tin,omegain: points where omega is given and the values there
64 c omega is assumed constant outside the given interval
65
66 open (20,file='status.log',form='formatted')
67 write (20,*) 1
68 close (20)
69
70 open (20,file='cvstag.log',form='formatted')
71 write (20,*) CVSTAG
72 close (20)
73
74 c read parameters etc.
75 call clnfile(5,10)
76 write (6,*) 'ierr, err'
77 read (10,*) ierr,err
78 write (6,*) 'name of file to read splittings from'
79 read (10,'(A80)') splitfile
80 write (6,*) 'control file for reading of eigenfunctions'
81 read (10,'(A80)') cefunc
82 write (6,*) 'number of points to skip in r for setting up f'
83 read (10,*) nskipr
84 write (6,*) 'number of points in x for setting up the plm''s'
85 read (10,*) nsetp
86 write (6,*) 'number of points to skip in x for setting up g'
87 read (10,*) nskipt
88 write (6,*) 'irdf, iwrf, irdg, iwrg, irdb, iwrb'
89 read (10,*) irdf,iwrf,irdg,iwrg,irdb,iwrb
90 write (6,*) '0 for multiple centerpoints, 1 for 1'
91 read (10,*) icenter
92 write (6,*) 'mur, mut'
93 read (10,*) mur,mut
94 nreg=2
95 write (6,*) 'icovar, icoeff, iavkern, nsr, nst, nsr2, nt2'
96 read (10,*) icovar, icoeff, iavkern, nsr, nst, nsr2, nt2
97 if (nsr.eq.0) then
98 nav=nst
99 do 105,i=1,nav
100 read (10,*) jav(i),kav(i)
101 105 continue
102 end if
103 c write (6,*) 'name of output file'
104 c read (10,'(A80)') outfile
105 close(10)
106
107 pi=4.0D0*atan(1.0D0)
108 pi2=pi/2.0D0
109 write (6,*) 'start reading splittings'
110 time1=dtime(dummy)
111 open (11,file=splitfile,status='old')
112 i=1
113 110 continue
114 if (ierr.eq.0) then
115 read (11,*,end=150) lnlm(i),nnlm1(i),fnlm(i),
116 c mnlm(i),itnlm(i),it1nlm(i),split(i)
117 sigma(i)=err
118 else
119 read (11,*,end=150) lnlm(i),nnlm1(i),fnlm(i),
120 c mnlm(i),itnlm(i),it1nlm(i),
121 c split(i),sigma(i)
122 end if
123 sigma1(i)=1.0/sigma(i)
124 ss1(i)=split(i)*sigma1(i)
125 ss2(i)=split(i)*sigma1(i)**2
126 i=i+1
127 goto 110
128 150 continue
129 close(11)
130 write (6,*) 'splittings have been read, time=',
131 c dtime(dummy),dummy(1),dummy(2)
132 nnlm=i-1
133 call setilm
134 write (6,*) dtime(dummy),dummy(1),dummy(2)
135 lmax1=-1
136 do 10,i=1,nnlm
137 if (lnlm(i).gt.lmax1) lmax1=lnlm(i)
138 10 continue
139 i1=0
140 do 100,l=0,lmax1
141 nl=0
142 do 20,i=1,nnlm
143 if (lnlm(i).eq.l) then
144 nl=nl+1
145 llist(nl)=i
146 end if
147 20 continue
148 i2=i1+1
149 do 60,i=1,nl
150 il=llist(i)
151 found=.false.
152 do 30,ii=i2,i1
153 if (nnlm1(il).eq.nmode(ii)) then
154 found=.true.
155 ifound=ii
156 end if
157 30 continue
158 if (found) then
159 iln(il)=ifound
160 else
161 i1=i1+1
162 lmode(i1)=l
163 nmode(i1)=nnlm1(il)
164 freq(i1)=fnlm(il)
165 iln(il)=i1
166 end if
167 60 continue
168 100 continue
169 nmodes=i1
170 write (6,*) dtime(dummy),dummy(1),dummy(2)
171 write (6,*) '#nl=',nmodes
172 write (6,*) '#nlm=',nnlm
173 c open (11,file=splitfile,form='unformatted',status='old')
174 c i=1
175 c i1=0
176 c110 read (11,end=150) lmode(i),nmode(i),freq(i),
177 c c (spliti(m),m=1,lmode(i))
178 c m0index(i)=i1
179 c do 130,m=1,lmode(i)
180 c i1=i1+1
181 c split(i1)=spliti(m)
182 c sigma(i1)=1.0E0
183 c sigma1(i1)=1.0D0/sigma(i1)
184 c ss1(i1)=split(i1)*sigma1(i1)
185 c lsplit(i1)=lmode(i)
186 c msplit(i1)=m
187 c isplit(i1)=i
188 c130 continue
189 c i=i+1
190 c goto 110
191 c150 nmodes=i-1
192 c nnlm=i1
193 c close (11)
194 write (6,*) 'start reading eigenfunctions'
195 time1=dtime(dummy)
196 call readefun(cefunc)
197 c do i=1,nrad
198 c write (6,*) dradmesh(i)
199 c end do
200 write (6,*) 'eigenfunctions have been read, time=',
201 c dtime(dummy),dummy(1),dummy(2)
202 if (mod(nrad-1,nskipr).ne.0) then
203 write (6,*) 'skip to get inversion mesh in radius does not'
204 write (6,*) 'divide into number of kernel meshpoints-1'
205 stop
206 end if
207 if (mod(nrad-1,nsr2).ne.0) then
208 write (6,*) 'skip to get av. kernel mesh in radius does not'
209 write (6,*) 'divide into number of kernel meshpoints-1'
210 stop
211 end if
212 nrad1=(nrad-1)/nskipr+1
213 nr2=(nrad-1)/nsr2+1
214 do 160,i=1,nrad1
215 rmesh1(i)=radmesh(nskipr*(i-1)+1)
216 160 continue
217 call dsetint(1,nrad1,rmesh1,rweights)
218 if (irdf.eq.0) then
219 write (6,*) 'start setting up f matrix'
220 time1=dtime(dummy)
221 call setf(nskipr,f1,f2)
222 write (6,*) 'f matrix has been set up, time=',
223 c dtime(dummy),dummy(1),dummy(2)
224 else
225 write (6,*) 'start reading f matrix'
226 time1=dtime(dummy)
227 open (21,file=tmpdir//'savef',form='unformatted',
228 c status='old')
229 read (21) ((f1(i,j),i=1,nmodes),j=1,nrad1)
230 read (21) ((f2(i,j),i=1,nmodes),j=1,nrad1)
231 close(21)
232 write (6,*) 'f matrix has been read, time=',
233 c dtime(dummy),dummy(1),dummy(2)
234 end if
235 if (iwrf.ne.0) then
236 write (6,*) 'start writing f matrix'
237 time1=dtime(dummy)
238 open (21,file=tmpdir//'savef',form='unformatted',
239 c status='unknown')
240 write (21) ((f1(i,j),i=1,nmodes),j=1,nrad1)
241 write (21) ((f2(i,j),i=1,nmodes),j=1,nrad1)
242 close(21)
243 write (6,*) 'f matrix has been written, time=',
244 c dtime(dummy),dummy(1),dummy(2)
245 end if
246 if (mod(nsetp-1,nskipt).ne.0) then
247 write (6,*) 'skip to get av. kernel mesh in latitude does not'
248 write (6,*) 'divide into original number of meshpoints-1'
249 stop
250 end if
251 nsetp1=(nsetp-1)/nskipt+1
252 if (irdg.eq.0) then
253 write (6,*) 'start setting up g matrix'
254 time1=dtime(dummy)
255 call setg(0,nsetp,nsetp1,g1,g2,tmesh1)
256 write (6,*) 'g matrix has been set up, time=',
257 c dtime(dummy),dummy(1),dummy(2)
258 else
259 write (6,*) 'start reading g matrix'
260 time1=dtime(dummy)
261 open (21,file=tmpdir//'saveg',form='unformatted',
262 c status='old')
263 read (21) (tmesh1(i),i=1,nsetp1)
264 c read (21) nlm,(lplm(i),mplm(i),i=1,nlm)
265 cc read (21) lmax,((iplm(i,j),i=0,lmax),j=0,lmax)
266 read (21) ((g1(i,j),i=1,nlm),j=1,nsetp1)
267 read (21) ((g2(i,j),i=1,nlm),j=1,nsetp1)
268 close(21)
269 write (6,*) 'g matrix has been read, time=',
270 c dtime(dummy),dummy(1),dummy(2)
271 end if
272 if (iwrg.ne.0) then
273 write (6,*) 'start writing g matrix'
274 time1=dtime(dummy)
275 open (21,file=tmpdir//'saveg',form='unformatted',
276 c status='unknown')
277 write (21) (tmesh1(i),i=1,nsetp1)
278 c write (21) nlm,(lplm(i),mplm(i),i=1,nlm)
279 cc write (21) lmax,((iplm(i,j),i=0,lmax),j=0,lmax)
280 write (21) ((g1(i,j),i=1,nlm),j=1,nsetp1)
281 write (21) ((g2(i,j),i=1,nlm),j=1,nsetp1)
282 close(21)
283 write (6,*) 'g matrix has been written, time=',
284 c dtime(dummy),dummy(1),dummy(2)
285 end if
286 call dsetint(1,nsetp1,tmesh1,tweights)
287 npoints=nrad1*nsetp1
288 nident=icenter*(nsetp1-1)
289 npoi1=npoints-nident
290 do 500,i2=1,npoints
291 k2=mod(i2-1,nsetp1)+1
292 j2=(i2-1)/nsetp1+1
293 ijk(j2,k2)=max(1,i2-nident)
294 500 continue
295 do 600,i=1,nnlm
296 c lmsplit(i)=iplm(msplit(i),lsplit(i))
297 600 continue
298 open (24,file='rmesh.orig',form='formatted')
299 do 610,i=1,nrad
300 write (24,*) radmesh(i)
301 610 continue
302 close (24)
303 if (iavkern.ne.0) then
304 open (24,file='avkern.log',form='formatted')
305 end if
306 open (25,file='rot.2d',form='formatted')
307 if (icovar.eq.1) then
308 open (26,file='err.2d',form='formatted')
309 end if
310 open (27,file='avkern.2d',form='unformatted')
311 if (icoeff.eq.1) then
312 open (28,file='coeff.2d',form='unformatted')
313 end if
314 xr(1)=-9.5
315 xr(2)=-8.5
316 xr(3)=-3.5
317 xr(4)=-3.0
318 xr(5)=-2.5
319 xr(6)=-2.0
320 c do 5000,imr=-8,-4,1
321 c mur=10.**imr
322 c do 5000,imr=1,6
323 c mur=10.**xr(imr)
324 c do 5000,imt=-7,-1,1
325 c mut=10.**imt
326 c write (6,*) mur,mut
327 if (irdb.eq.0) then
328 write (6,*) 'start setting up b matrix'
329 time1=dtime(dummy)
330 do 950,i=1,npoi1
331 do 940,j=1,npoi1
332 b(j,i)=0.0D0
333 940 continue
334 950 continue
335 nb1=0
336 nb2=0
337 nb2a=0
338 t2=0
339 t2a=0
340 do 1100,is1=1,nnlm,maxs1
341 write (6,*) is1
342 nbig=0
343 nbig2=0
344 is0=min(maxs1-1,nnlm-is1)
345 ix=0
346 do 1020,i2=1,npoints
347 ibig(i2)=.false.
348 k2=mod(i2-1,nsetp1)+1
349 j2=(i2-1)/nsetp1+1
350 nbig1=0
351 do 1010,isx=0,is0
352 is=is1+isx
353 i=iln(is)
354 ilmi=ilm(is)
355 ax(isx)=sigma1(is)*
356 c (f1(i,j2)*g1(ilmi,k2)+f2(i,j2)*g2(ilmi,k2))
357 if (abs(ax(isx)).gt.1e-15) then
358 nbig1=nbig1+1
359 end if
360 1010 continue
361 if (nbig1.ne.0) then
362 ibig(i2)=.true.
363 nbig=nbig+nbig1
364 nbig2=nbig2+1
365 ix=ix+1
366 iax(ix)=i2
367 do 1015,isx=0,is0
368 a(isx,ix)=ax(isx)
369 1015 continue
370 end if
371 1020 continue
372 nx=ix
373 nb1=nb1+nbig
374 nb2=nb2+nbig2
375 nb2a=nb2a+npoints
376 t2=t2+nbig2**2*(is0+1.)
377 t2a=t2a+npoints**2*(is0+1.)
378 do 1080,i1x=1,nx,nblock
379 n1x=min(i1x+nblock-1,nx)-i1x+1
380 do 1070,i2x=i1x,nx,nblock
381 n2x=min(i2x+nblock-1,nx)-i2x+1
382 call dgemm('t','n',n2x,n1x,is0+1,1.0D0,a(0,i2x),maxs1,
383 c a(0,i1x),maxs1,0.0D0,bx,nblock)
384 do 1060,i1=i1x,min(i1x+nblock-1,nx)
385 i1a=max(1,iax(i1)-nident)
386 do 1050,i2=i2x,min(i2x+nblock-1,nx)
387 i2a=max(1,iax(i2)-nident)
388 if (i2a.ge.i1a) then
389 c b(i2a,i1a)=b(i2a,i1a)+ddot(is0+1,a(0,i2),1,a(0,i1),1)
390 b(i2a,i1a)=b(i2a,i1a)+bx(i2-i2x+1,i1-i1x+1)
391 end if
392 1050 continue
393 1060 continue
394 1070 continue
395 1080 continue
396 1100 continue
397 write (6,*) (nb1+0.0)/(nnlm*npoints),(nb2+0.0)/nb2a,t2/t2a
398 do 1290,i2=1,npoi1
399 do 1280,i1=1,i2-1
400 b(i1,i2)=b(i2,i1)
401 1280 continue
402 1290 continue
403 write (6,*) 'b matrix has been set up, time=',
404 c dtime(dummy),dummy(1),dummy(2)
405 else
406 write (6,*) 'start reading b matrix'
407 time1=dtime(dummy)
408 open (21,file=tmpdir//'saveb',form='unformatted',
409 c status='old')
410 read (21) ((b(j,i),j=1,npoi1),i=1,npoi1)
411 close(21)
412 c do i=1,npoi1
413 c do j=1,npoi1
414 c b(j,i)=real(b(j,i))
415 c end do
416 c end do
417 write (6,*) 'b matrix has been read, time=',
418 c dtime(dummy),dummy(1),dummy(2)
419 end if
420 c write (6,*) (b(i,i),i=1,npoi1)
421 if (iwrb.ne.0) then
422 write (6,*) 'start writing b matrix'
423 time1=dtime(dummy)
424 open (21,file=tmpdir//'saveb',form='unformatted',
425 c status='unknown')
426 write (21) ((b(j,i),j=1,npoi1),i=1,npoi1)
427 close(21)
428 write (6,*) 'b matrix has been written, time=',
429 c dtime(dummy),dummy(1),dummy(2)
430 c to prevent recalculating b if doing trade-off curves
431 irdb=1
432 end if
433 write (6,*) 'start adding regularization'
434 time1=dtime(dummy)
435 do 1310,i=1,npoi1
436 do 1300,j=1,npoi1
437 d(j,i)=b(j,i)
438 1300 continue
439 1310 continue
440 c c1 is is the regularization in the 1-d radial case
441 c times the radius
442 if (nreg.eq.2) then
443 c set up weighting functions
444 do 1400,i=1,nrad1
445 c fr(i)=1.0
446 fr(i)=rmesh1(i)
447 if (rmesh1(i).eq.0) then
448 ft(i)=0.0
449 else
450 ft(i)=1.0D0/rmesh1(i)**4
451 end if
452 c just a try
453 c ft(i)=1.0
454 1400 continue
455 do 1405,i=1,nrad1
456 fwr(i)=mur*rweights(i)*fr(i)
457 c fwr(i)=mur*rweights(i)*rmesh1(i)
458 1405 continue
459 do 1430,j=1,nrad1
460 do 1420,k=j,min(nrad1,j+2)
461 ckj=0.0D0
462 do 1410,i=max(2,j-1,k-1),min(nrad1-1,j+1,k+1)
463 xm1=rmesh1(i-1)-rmesh1(i)
464 xp1=rmesh1(i+1)-rmesh1(i)
465 den=xm1*xp1*(xm1-xp1)
466 if (i.eq.j-1) fji=-xm1
467 if (i.eq.j) fji=xm1-xp1
468 if (i.eq.j+1) fji=xp1
469 if (i.eq.k-1) fki=-xm1
470 if (i.eq.k) fki=xm1-xp1
471 if (i.eq.k+1) fki=xp1
472 ckj=ckj+fji*fki*fwr(i)/den**2
473 1410 continue
474 c factor of 4=2^2 since second derivative is 2*a where y=a x^2+b x +c
475 c1(k,j)=4.*ckj
476 c1(j,k)=4.*ckj
477 1420 continue
478 1430 continue
479 do 1500,i=1,nsetp1
480 fwt(i)=mut*tweights(i)
481 1500 continue
482 do 1530,j=1,nsetp1
483 do 1520,k=j,min(nsetp1,j+2)
484 ckj=0.0D0
485 do 1510,i=max(1,j-1,k-1),min(nsetp1,j+1,k+1)
486 t=tmesh1(i)
487 if (i.gt.1) tm1=tmesh1(i-1)
488 if (i.eq.1) tm1=2*tmesh1(1)-tmesh1(2)
489 if (i.lt.nsetp1) tp1=tmesh1(i+1)
490 if (i.eq.nsetp1) tp1=2*tmesh1(nsetp1)-tmesh1(nsetp1-1)
491 xm1=tm1-t
492 xp1=tp1-t
493 den=xm1*xp1*(xm1-xp1)
494 if (i.eq.j-1) fji=-xm1
495 if (i.eq.j) fji=xm1-xp1
496 if (i.eq.j+1) fji=xp1
497 if (i.eq.1.and.j.eq.2) fji=xp1-xm1
498 if (i.eq.nsetp1.and.j.eq.(nsetp1-1)) fji=xp1-xm1
499 if (i.eq.k-1) fki=-xm1
500 if (i.eq.k) fki=xm1-xp1
501 if (i.eq.k+1) fki=xp1
502 if (i.eq.1.and.k.eq.2) fki=xp1-xm1
503 if (i.eq.nsetp1.and.k.eq.(nsetp1-1)) fki=xp1-xm1
504 ckj=ckj+fji*fki*fwt(i)/den**2
505 c write (6,*) i,j,k,ckj
506 1510 continue
507 c factor of 4=2^2 since second derivative is 2*a where y=a x^2+b x +c
508 c2(k,j)=4.*ckj
509 c2(j,k)=4.*ckj
510 1520 continue
511 1530 continue
512 end if
513 write (6,*) dtime(dummy),dummy(1),dummy(2)
514 c first the regularization in theta is applied
515 do 1650,j=2,nrad1
516 wr3i=rweights(j)*ft(j)
517 c wr3i=rweights(j)/rmesh1(j)**4
518 do 1640,k1=1,nsetp1
519 i1=ijk(j,k1)
520 do 1630,k2=max(1,k1-2),min(nsetp1,k1+2)
521 i2=ijk(j,k2)
522 d(i2,i1)=d(i2,i1)+c2(k2,k1)*wr3i
523 1630 continue
524 1640 continue
525 1650 continue
526 write (6,*) dtime(dummy),dummy(1),dummy(2)
527 c then the regularization in radius is added
528 do 1750,k=1,nsetp1
529 w=tweights(k)
530 do 1740,j1=1,nrad1
531 i1=ijk(j1,k)
532 do 1730,j2=1,nrad1
533 i2=ijk(j2,k)
534 d(i2,i1)=d(i2,i1)+c1(j2,j1)*w
535 1730 continue
536 1740 continue
537 1750 continue
538 write (6,*) 'regularization has been added, time=',
539 c dtime(dummy),dummy(1),dummy(2)
540 write (6,*) 'start setting up right hand side'
541 time1=dtime(dummy)
542 do 1810,i=1,npoi1
543 right(i)=0.0D0
544 1810 continue
545 do 1850,i1=1,npoints
546 c i1a=max(1,i1-nident)
547 k1=mod(i1-1,nsetp1)+1
548 j1=(i1-1)/nsetp1+1
549 do 1820,i=1,nmodes
550 f1x(i)=f1(i,j1)
551 f2x(i)=f2(i,j1)
552 1820 continue
553 do 1830,i=1,nlm
554 g1x(i)=g1(i,k1)
555 g2x(i)=g2(i,k1)
556 1830 continue
557 r=0.0D0
558 do 1840,is=1,nnlm
559 i=iln(is)
560 ilmi=ilm(is)
561 r=r+ss2(is)*(f1x(i)*g1x(ilmi)+f2x(i)*g2x(ilmi))
562 1840 continue
563 c right(i1a)=right(i1a)+r
564 help(i1)=r
565 1850 continue
566 do 1860,i1=1,npoints
567 i1a=max(1,i1-nident)
568 right(i1a)=right(i1a)+help(i1)
569 1860 continue
570 write (6,*) 'right hand side has been set up, time=',
571 c dtime(dummy),dummy(1),dummy(2)
572 c open (21,file=tmpdir//'saved',form='unformatted',
573 c c status='unknown')
574 c do i=1,npoi1
575 c write (21) (d(j,i),j=1,npoi1)
576 c end do
577 c write (21) (right(i),i=1,npoi1)
578 c close(21)
579 c stop
580 write (6,*) 'start decomposing d matrix'
581 time1=dtime(dummy)
582 c call dpoco(d,maxpoints,npoi1,rcond,help,info)
583 anorm=dlansy('1','U',npoi1,d,maxpoints,help)
584 call dpotrf('U',npoi1,d,maxpoints,info)
585 write (6,*) 'd matrix has been decomposed, time=',
586 c dtime(dummy),dummy(1),dummy(2),info
587 write (6,*) 'start finding condition number'
588 call dpocon('U',npoi1,d,maxpoints,anorm,rcond,help,iwork,info1)
589 write (6,*) 'condition number has been found, time=',
590 c dtime(dummy),dummy(1),dummy(2)
591 write (6,*) rcond,info,info1,anorm
592 write (6,*) 'start finding solution'
593 time1=dtime(dummy)
594 c call dposl(d,maxpoints,npoi1,right)
595 call dpotrs('U',npoi1,1,d,maxpoints,right,maxpoints,info)
596 write (6,*) 'solution has been found, time=',
597 c dtime(dummy),dummy(1),dummy(2)
598 c write (21) (right(i),i=1,npoi1)
599 c close(21)
600 c stop
601 write (6,*) 'rcond,info',rcond,info
602 do 1875,j=1,nrad1
603 do 1870,k=1,nsetp1
604 omega(j,k)=right(ijk(j,k))
605 1870 continue
606 1875 continue
607 do 1877,j=1,nrad1
608 write (25,'(101f12.4)') (right(ijk(j,k)),k=1,nsetp1)
609 1877 continue
610 write (6,*) 'start calculating splittings'
611 time1=dtime(dummy)
612 c set ixl to the first mode with a given (l,n)
613 c set ixh to the last mode with a given (l,n)
614 ilast=-1
615 ix=0
616 do 1880,i1=1,nnlm
617 i=iln(i1)
618 if (i.ne.ilast) then
619 ix=ix+1
620 ixl(ix)=i1
621 ilast=i
622 end if
623 ixh(ix)=i1
624 1880 continue
625 nx=ix
626 chisq=0.0
627 c*$* assert do (serial)
628 open(10,file='splittings.out',form='formatted')
629 do 1897,ix=1,nx
630 c loop first over each ln
631 i=iln(ixl(ix))
632 do 1891,j=1,nrad1
633 f1i(j)=f1(i,j)
634 f2i(j)=f2(i,j)
635 1891 continue
636 call dgemv('t',nrad1,nsetp1,1.0D0,omega,maxr1,
637 c f1i,1,0.0D0,f1o,1)
638 call dgemv('t',nrad1,nsetp1,1.0D0,omega,maxr1,
639 c f2i,1,0.0D0,f2o,1)
640 c now do the individual nlm's
641 do 1895,i1=ixl(ix),ixh(ix)
642 index=ilm(i1)
643 c do 1893,k=1,nsetp1
644 c g1i(k)=g1(index,k)
645 c g2i(k)=g2(index,k)
646 c1893 continue
647 c s=ddot(nsetp1,f1o,1,g1i,1)+ddot(nsetp1,f2o,1,g2i,1)
648 s=0.0D0
649 do 1893,k=1,nsetp1
650 s=s+g1(index,k)*f1o(k)+g2(index,k)*f2o(k)
651 1893 continue
652 l=lmode(i)
653 n=nmode(i)
654 write (10,'(i4,i4,f8.2,i4,2i3,3f12.4)')
655 c write (6,'(i4,i4,f8.2,i4,2i3,3f12.4)')
656 c l,n,freq(i),mnlm(i1),itnlm(i1),it1nlm(i1),
657 c s,split(i1),sigma(i1)
658 chisq=chisq+((s-split(i1))*sigma1(i1))**2
659 1895 continue
660 1897 continue
661 write (6,*) 'splittings have been calculated, time=',
662 c dtime(dummy),dummy(1),dummy(2)
663 write (6,*) 'start calculating wiggle measures'
664 time1=dtime(dummy)
665 do 1907,k=1,nsetp1
666 do 1905,j=1,nrad1
667 fro(j,k)=0.0
668 fto(j,k)=0.0
669 c omega(j,k)=sin(tmesh1(k))**2*rmesh1(j)**2
670 c omega(j,k)=rmesh1(j)**3
671 1905 continue
672 1907 continue
673 do 1930,j=1,nrad1
674 do 1920,i=max(2,j-1),min(nrad1-1,j+1)
675 xm1=rmesh1(i-1)-rmesh1(i)
676 xp1=rmesh1(i+1)-rmesh1(i)
677 den=xm1*xp1*(xm1-xp1)
678 if (i.eq.j-1) fji=-xm1
679 if (i.eq.j) fji=xm1-xp1
680 if (i.eq.j+1) fji=xp1
681 do 1910,k=1,nsetp1
682 fro(i,k)=fro(i,k)+omega(j,k)*2*fji/den
683 1910 continue
684 1920 continue
685 1930 continue
686 do 1960,k=1,nsetp1
687 do 1950,i=max(1,k-1),min(nsetp1,k+1)
688 t=tmesh1(i)
689 if (i.gt.1) tm1=tmesh1(i-1)
690 if (i.eq.1) tm1=2*tmesh1(1)-tmesh1(2)
691 if (i.lt.nsetp1) tp1=tmesh1(i+1)
692 if (i.eq.nsetp1) tp1=2*tmesh1(nsetp1)-tmesh1(nsetp1-1)
693 xm1=tm1-t
694 xp1=tp1-t
695 den=xm1*xp1*(xm1-xp1)
696 if (i.eq.k-1) fki=-xm1
697 if (i.eq.k) fki=xm1-xp1
698 if (i.eq.k+1) fki=xp1
699 if (i.eq.1.and.k.eq.2) fki=xp1-xm1
700 if (i.eq.nsetp1.and.k.eq.(nsetp1-1)) fki=xp1-xm1
701 do 1940,j=1,nrad1
702 fto(j,i)=fto(j,i)+omega(j,k)*2*fki/den
703 1940 continue
704 1950 continue
705 1960 continue
706 c calculate wigr and wigt
707 wigr=0.0
708 wigt=0.0
709 do 1980,k=1,nsetp1
710 wt=tweights(k)
711 do 1970,j=1,nrad1
712 wigr=wigr+wt*rweights(j)*fr(j)*fro(j,k)**2
713 wigt=wigt+wt*rweights(j)*ft(j)*fto(j,k)**2
714 1970 continue
715 1980 continue
716 write (6,*) 'wiggle measures have been computed, time=',
717 c dtime(dummy),dummy(1),dummy(2)
718 write (6,*) 'chisq=',chisq
719 write (6,*) 'rms=',sqrt(chisq/nnlm)
720 write (6,*) 'wigr, wigt=',wigr,wigt
721 write (6,*) '%',mur,mut,chisq,sqrt(chisq/nnlm),wigr,wigt,
722 c info,rcond
723 c end do
724 c end do
725 c stop
726 if ((icovar.eq.1).or.(iavkern.ge.1)) then
727 write (6,*) 'start inverting d matrix'
728 time1=dtime(dummy)
729 c call dpodi(d,maxpoints,npoi1,det,1)
730 call dpotri('U',npoi1,d,maxpoints,info)
731 c open (21,file=tmpdir//'saved1',form='unformatted',
732 c c status='unknown',access='direct',recl=8*npoi1)
733 do 2095,i=1,npoi1
734 c do 2092,j=1,npoi1
735 c di(j)=0.0D0
736 c2092 continue
737 c di(i)=1.0D0
738 c call dposl(d,maxpoints,npoi1,di)
739 c write (21,rec=i) (di(j),j=1,npoi1)
740 c write (21,rec=i) (d(j,i),j=1,i),(d(i,j),j=i+1,npoi1)
741 do 2092,j=i+1,npoi1
742 d(j,i)=d(i,j)
743 2092 continue
744 2095 continue
745 c close(21)
746 write (6,*) 'd matrix has been inverted, time=',
747 c dtime(dummy),dummy(1),dummy(2)
748 end if
749 if (icovar.ge.1) then
750 write (6,*) 'start finding variances'
751 c the following loop could be blocked
752 c do 2150,i=1,npoi1
753 c call dgemv('n',npoi1,npoi1,1.0D0,b,maxpoints,d(1,i),1,
754 c c 0.0D0,dib,1)
755 c var(i)=ddot(npoi1,d(1,i),1,dib,1)
756 c2150 continue
757 do 2150,i=1,npoi1,nblock1
758 ix=min(npoi1,i+nblock1-1)
759 nx=ix+1-i
760 call dgemm('n','n',npoi1,nx,npoi1,1.0D0,b,maxpoints,
761 c d(1,i),maxpoints,0.0D0,dibx,maxpoints)
762 do 2140,j=i,ix
763 var(j)=ddot1(npoi1,d(1,j),dibx(1,j-i+1))
764 2140 continue
765 2150 continue
766 do 2190,j=1,nrad1
767 write (26,'(100f11.5)') (sqrt(var(ijk(j,k))),k=1,nsetp1)
768 2190 continue
769 c close(21)
770 write (6,*) 'variances have been found, time=',
771 c dtime(dummy),dummy(1),dummy(2)
772 c else
773 c write (6,*) 'start finding covariance matrix'
774 c time1=dtime(dummy)
775 c do 2230,i=1,npoi1
776 c do 2220,j=1,npoi1
777 c d(j,i)=ddot1(npoi1,b(1,j),d1(1,i))
778 c2220 continue
779 c2230 continue
780 c do 2250,i=1,npoi1
781 c do 2240,j=i,npoi1
782 c cov(j,i)=ddot1(npoi1,d1(1,j),d(1,i))
783 c2240 continue
784 c2250 continue
785 c do 2270,i=1,npoi1
786 c do 2260,j=1,i-1
787 c cov(j,i)=cov(i,j)
788 c2260 continue
789 c2270 continue
790 c do 2280,i=1,nrad1
791 c write (6,*) rmesh1(i)
792 c2280 continue
793 c do 2290,j=1,nrad1
794 c write (6,'(100f9.5)') (sqrt(cov(ijk(j,k),ijk(j,k))),
795 c c k=1,nsetp1)
796 c2290 continue
797 c write (6,*) 'covariance matrix has been found, time=',
798 c c dtime(dummy),dummy(1),dummy(2)
799 end if
800 if (iavkern.eq.0) then
801 write (6,*) 'skipping calculation of averaging kernels'
802 else
803 write (6,*) 'start finding averaging kernels'
804 time1=dtime(dummy)
805 write (6,*) nt2,nr2
806 c open (21,file=tmpdir//'saved1',form='unformatted',
807 c c status='old',access='direct',recl=8*npoi1)
808 c moved up
809 c open (27,file='avkern.2d',form='unformatted')
810 c open (28,file='coeff',form='unformatted')
811 do 2300,j1=1,nr2
812 x2(j1)=radmesh(nsr2*(j1-1)+1)
813 2300 continue
814 pi=4.0D0*atan(1.0D0)
815 pi2=pi/2.0D0
816 do 2301,k1=1,nt2
817 x1(k1)=pi2*(k1-1)/(nt2-1.)
818 2301 end do
819 write (6,*) nrad1,nsr,nsetp1,nst
820 if (nsr.ne.0) then
821 i=0
822 do 2303,j=1,nrad1,nsr
823 do 2302,k=1,nsetp1,nst
824 i=i+1
825 jav(i)=j
826 kav(i)=k
827 2302 continue
828 2303 continue
829 nav=i
830 end if
831 c do 2380,j=1,nrad1,nsr
832 c do 2370,k=1,nsetp1,nst
833 do 2380,iav=1,nav
834 j=jav(iav)
835 k=kav(iav)
836 ipoi=ijk(j,k)
837 c read (21,rec=ipoi) (di(j1),j1=1,npoi1)
838 do 2305,is=1,nnlm
839 coeff(is)=0.0D0
840 2305 continue
841 do 2330,i1=1,npoints
842 i1a=max(1,i1-nident)
843 k1=mod(i1-1,nsetp1)+1
844 j1=(i1-1)/nsetp1+1
845 do 2310,i=1,nmodes
846 f1x(i)=f1(i,j1)
847 f2x(i)=f2(i,j1)
848 2310 continue
849 do 2315,i=1,nlm
850 g1x(i)=g1(i,k1)
851 g2x(i)=g2(i,k1)
852 2315 continue
853 c d1i=d1(i1a,ipoi)
854 c d1i=di(i1a)
855 d1i=d(i1a,ipoi)
856 do 2320,is=1,nnlm
857 i=iln(is)
858 ilmi=ilm(is)
859 coeff(is)=coeff(is)+d1i*sigma1(is)*
860 c (f1x(i)*g1x(ilmi)+f2x(i)*g2x(ilmi))
861 2320 continue
862 2330 continue
863 c s1=0.0D0
864 c s2=0.0D0
865 c s3=0.0D0
866 vari=0.0
867 do 2350,is=1,nnlm
868 coeff(is)=coeff(is)*sigma1(is)
869 c s1=s1+coeff(is)*split(is)
870 c s2=s2+mnlm(is)*coeff(is)
871 c s3=s3+coeff(is)
872 vari=vari+(coeff(is)*sigma(is))**2
873 2350 continue
874 c write (6,*) '#',j,k
875 c write (6,*) s1,right(ipoi),s1-right(ipoi),s2,s3
876 c call set2davk(nnlm,isplit,msplit,coeff,avkern,nsr2,nt2)
877 call set2davk(coeff,avkern,nsr2,nt2)
878 if (iavkern.eq.2) then
879 write (27) ((avkern(k1,j1),k1=1,nt2),j1=1,nr2)
880 end if
881 c call fitmn(x1,x2,nt2,nr2,avkern,yy2,mjta,4,
882 c c real(tmesh1(k)),real(rmesh1(j)))
883 write (24,*) mur,mut,real(rmesh1(j)),real(tmesh1(k)),
884 c right(ipoi),sqrt(vari),
885 c mjta(3),mjta(1),mjta(4),mjta(2)
886 if (icoeff.eq.1) then
887 write (28) (real(coeff(is)),is=1,nnlm)
888 end if
889 c end if
890 c2370 continue
891 2380 continue
892 c close(27)
893 c close(28)
894 c close(21)
895 write (6,*) 'averaging kernels set up, time=',
896 c dtime(dummy),dummy(1),dummy(2)
897 end if
898 5000 continue
899 if (iavkern.ne.0) then
900 close (24)
901 end if
902 close (25)
903 if (icovar.eq.26) then
904 close (26)
905 end if
906 close (27)
907 if (icoeff.eq.28) then
908 close (28)
909 end if
910 open (20,file='status.log',form='formatted')
911 write (20,*) 0
912 close (20)
913
914 write (6,*) 'total time used:',etime(dummy)
915 write (6,*) 'user time: ',dummy(1)
916 write (6,*) 'system time: ',dummy(2)
917 end
918