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