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 if (icoeff.eq.1) then
307 open (28,file='coeff.2d',form='unformatted')
308 end if
309 xr(1)=-9.5
310 xr(2)=-8.5
311 xr(3)=-3.5
312 xr(4)=-3.0
313 xr(5)=-2.5
314 xr(6)=-2.0
315 c do 5000,imr=-8,-4,1
316 tplarson 1.1 c mur=10.**imr
317 c do 5000,imr=1,6
318 c mur=10.**xr(imr)
319 c do 5000,imt=-7,-1,1
320 c mut=10.**imt
321 c write (6,*) mur,mut
322 if (irdb.eq.0) then
323 write (6,*) 'start setting up b matrix'
324 time1=dtime(dummy)
325 do 950,i=1,npoi1
326 do 940,j=1,npoi1
327 b(j,i)=0.0D0
328 940 continue
329 950 continue
330 nb1=0
331 nb2=0
332 nb2a=0
333 t2=0
334 t2a=0
335 do 1100,is1=1,nnlm,maxs1
336 write (6,*) is1
337 tplarson 1.1 nbig=0
338 nbig2=0
339 is0=min(maxs1-1,nnlm-is1)
340 ix=0
341 do 1020,i2=1,npoints
342 ibig(i2)=.false.
343 k2=mod(i2-1,nsetp1)+1
344 j2=(i2-1)/nsetp1+1
345 nbig1=0
346 do 1010,isx=0,is0
347 is=is1+isx
348 i=iln(is)
349 ilmi=ilm(is)
350 ax(isx)=sigma1(is)*
351 c (f1(i,j2)*g1(ilmi,k2)+f2(i,j2)*g2(ilmi,k2))
352 if (abs(ax(isx)).gt.1e-15) then
353 nbig1=nbig1+1
354 end if
355 1010 continue
356 if (nbig1.ne.0) then
357 ibig(i2)=.true.
358 tplarson 1.1 nbig=nbig+nbig1
359 nbig2=nbig2+1
360 ix=ix+1
361 iax(ix)=i2
362 do 1015,isx=0,is0
363 a(isx,ix)=ax(isx)
364 1015 continue
365 end if
366 1020 continue
367 nx=ix
368 nb1=nb1+nbig
369 nb2=nb2+nbig2
370 nb2a=nb2a+npoints
371 t2=t2+nbig2**2*(is0+1.)
372 t2a=t2a+npoints**2*(is0+1.)
373 do 1080,i1x=1,nx,nblock
374 n1x=min(i1x+nblock-1,nx)-i1x+1
375 do 1070,i2x=i1x,nx,nblock
376 n2x=min(i2x+nblock-1,nx)-i2x+1
377 call dgemm('t','n',n2x,n1x,is0+1,1.0D0,a(0,i2x),maxs1,
378 c a(0,i1x),maxs1,0.0D0,bx,nblock)
379 tplarson 1.1 do 1060,i1=i1x,min(i1x+nblock-1,nx)
380 i1a=max(1,iax(i1)-nident)
381 do 1050,i2=i2x,min(i2x+nblock-1,nx)
382 i2a=max(1,iax(i2)-nident)
383 if (i2a.ge.i1a) then
384 c b(i2a,i1a)=b(i2a,i1a)+ddot(is0+1,a(0,i2),1,a(0,i1),1)
385 b(i2a,i1a)=b(i2a,i1a)+bx(i2-i2x+1,i1-i1x+1)
386 end if
387 1050 continue
388 1060 continue
389 1070 continue
390 1080 continue
391 1100 continue
392 write (6,*) (nb1+0.0)/(nnlm*npoints),(nb2+0.0)/nb2a,t2/t2a
393 do 1290,i2=1,npoi1
394 do 1280,i1=1,i2-1
395 b(i1,i2)=b(i2,i1)
396 1280 continue
397 1290 continue
398 write (6,*) 'b matrix has been set up, time=',
399 c dtime(dummy),dummy(1),dummy(2)
400 tplarson 1.1 else
401 write (6,*) 'start reading b matrix'
402 time1=dtime(dummy)
403 open (21,file=tmpdir//'saveb',form='unformatted',
404 c status='old')
405 read (21) ((b(j,i),j=1,npoi1),i=1,npoi1)
406 close(21)
407 c do i=1,npoi1
408 c do j=1,npoi1
409 c b(j,i)=real(b(j,i))
410 c end do
411 c end do
412 write (6,*) 'b matrix has been read, time=',
413 c dtime(dummy),dummy(1),dummy(2)
414 end if
415 c write (6,*) (b(i,i),i=1,npoi1)
416 if (iwrb.ne.0) then
417 write (6,*) 'start writing b matrix'
418 time1=dtime(dummy)
419 open (21,file=tmpdir//'saveb',form='unformatted',
420 c status='unknown')
421 tplarson 1.1 write (21) ((b(j,i),j=1,npoi1),i=1,npoi1)
422 close(21)
423 write (6,*) 'b matrix has been written, time=',
424 c dtime(dummy),dummy(1),dummy(2)
425 c to prevent recalculating b if doing trade-off curves
426 irdb=1
427 end if
428 write (6,*) 'start adding regularization'
429 time1=dtime(dummy)
430 do 1310,i=1,npoi1
431 do 1300,j=1,npoi1
432 d(j,i)=b(j,i)
433 1300 continue
434 1310 continue
435 c c1 is is the regularization in the 1-d radial case
436 c times the radius
437 if (nreg.eq.2) then
438 c set up weighting functions
439 do 1400,i=1,nrad1
440 c fr(i)=1.0
441 fr(i)=rmesh1(i)
442 tplarson 1.1 if (rmesh1(i).eq.0) then
443 ft(i)=0.0
444 else
445 ft(i)=1.0D0/rmesh1(i)**4
446 end if
447 c just a try
448 c ft(i)=1.0
449 1400 continue
450 do 1405,i=1,nrad1
451 fwr(i)=mur*rweights(i)*fr(i)
452 c fwr(i)=mur*rweights(i)*rmesh1(i)
453 1405 continue
454 do 1430,j=1,nrad1
455 do 1420,k=j,min(nrad1,j+2)
456 ckj=0.0D0
457 do 1410,i=max(2,j-1,k-1),min(nrad1-1,j+1,k+1)
458 xm1=rmesh1(i-1)-rmesh1(i)
459 xp1=rmesh1(i+1)-rmesh1(i)
460 den=xm1*xp1*(xm1-xp1)
461 if (i.eq.j-1) fji=-xm1
462 if (i.eq.j) fji=xm1-xp1
463 tplarson 1.1 if (i.eq.j+1) fji=xp1
464 if (i.eq.k-1) fki=-xm1
465 if (i.eq.k) fki=xm1-xp1
466 if (i.eq.k+1) fki=xp1
467 ckj=ckj+fji*fki*fwr(i)/den**2
468 1410 continue
469 c factor of 4=2^2 since second derivative is 2*a where y=a x^2+b x +c
470 c1(k,j)=4.*ckj
471 c1(j,k)=4.*ckj
472 1420 continue
473 1430 continue
474 do 1500,i=1,nsetp1
475 fwt(i)=mut*tweights(i)
476 1500 continue
477 do 1530,j=1,nsetp1
478 do 1520,k=j,min(nsetp1,j+2)
479 ckj=0.0D0
480 do 1510,i=max(1,j-1,k-1),min(nsetp1,j+1,k+1)
481 t=tmesh1(i)
482 if (i.gt.1) tm1=tmesh1(i-1)
483 if (i.eq.1) tm1=2*tmesh1(1)-tmesh1(2)
484 tplarson 1.1 if (i.lt.nsetp1) tp1=tmesh1(i+1)
485 if (i.eq.nsetp1) tp1=2*tmesh1(nsetp1)-tmesh1(nsetp1-1)
486 xm1=tm1-t
487 xp1=tp1-t
488 den=xm1*xp1*(xm1-xp1)
489 if (i.eq.j-1) fji=-xm1
490 if (i.eq.j) fji=xm1-xp1
491 if (i.eq.j+1) fji=xp1
492 if (i.eq.1.and.j.eq.2) fji=xp1-xm1
493 if (i.eq.nsetp1.and.j.eq.(nsetp1-1)) fji=xp1-xm1
494 if (i.eq.k-1) fki=-xm1
495 if (i.eq.k) fki=xm1-xp1
496 if (i.eq.k+1) fki=xp1
497 if (i.eq.1.and.k.eq.2) fki=xp1-xm1
498 if (i.eq.nsetp1.and.k.eq.(nsetp1-1)) fki=xp1-xm1
499 ckj=ckj+fji*fki*fwt(i)/den**2
500 c write (6,*) i,j,k,ckj
501 1510 continue
502 c factor of 4=2^2 since second derivative is 2*a where y=a x^2+b x +c
503 c2(k,j)=4.*ckj
504 c2(j,k)=4.*ckj
505 tplarson 1.1 1520 continue
506 1530 continue
507 end if
508 write (6,*) dtime(dummy),dummy(1),dummy(2)
509 c first the regularization in theta is applied
510 do 1650,j=2,nrad1
511 wr3i=rweights(j)*ft(j)
512 c wr3i=rweights(j)/rmesh1(j)**4
513 do 1640,k1=1,nsetp1
514 i1=ijk(j,k1)
515 do 1630,k2=max(1,k1-2),min(nsetp1,k1+2)
516 i2=ijk(j,k2)
517 d(i2,i1)=d(i2,i1)+c2(k2,k1)*wr3i
518 1630 continue
519 1640 continue
520 1650 continue
521 write (6,*) dtime(dummy),dummy(1),dummy(2)
522 c then the regularization in radius is added
523 do 1750,k=1,nsetp1
524 w=tweights(k)
525 do 1740,j1=1,nrad1
526 tplarson 1.1 i1=ijk(j1,k)
527 do 1730,j2=1,nrad1
528 i2=ijk(j2,k)
529 d(i2,i1)=d(i2,i1)+c1(j2,j1)*w
530 1730 continue
531 1740 continue
532 1750 continue
533 write (6,*) 'regularization has been added, time=',
534 c dtime(dummy),dummy(1),dummy(2)
535 write (6,*) 'start setting up right hand side'
536 time1=dtime(dummy)
537 do 1810,i=1,npoi1
538 right(i)=0.0D0
539 1810 continue
540 do 1850,i1=1,npoints
541 c i1a=max(1,i1-nident)
542 k1=mod(i1-1,nsetp1)+1
543 j1=(i1-1)/nsetp1+1
544 do 1820,i=1,nmodes
545 f1x(i)=f1(i,j1)
546 f2x(i)=f2(i,j1)
547 tplarson 1.1 1820 continue
548 do 1830,i=1,nlm
549 g1x(i)=g1(i,k1)
550 g2x(i)=g2(i,k1)
551 1830 continue
552 r=0.0D0
553 do 1840,is=1,nnlm
554 i=iln(is)
555 ilmi=ilm(is)
556 r=r+ss2(is)*(f1x(i)*g1x(ilmi)+f2x(i)*g2x(ilmi))
557 1840 continue
558 c right(i1a)=right(i1a)+r
559 help(i1)=r
560 1850 continue
561 do 1860,i1=1,npoints
562 i1a=max(1,i1-nident)
563 right(i1a)=right(i1a)+help(i1)
564 1860 continue
565 write (6,*) 'right hand side has been set up, time=',
566 c dtime(dummy),dummy(1),dummy(2)
567 c open (21,file=tmpdir//'saved',form='unformatted',
568 tplarson 1.1 c c status='unknown')
569 c do i=1,npoi1
570 c write (21) (d(j,i),j=1,npoi1)
571 c end do
572 c write (21) (right(i),i=1,npoi1)
573 c close(21)
574 c stop
575 write (6,*) 'start decomposing d matrix'
576 time1=dtime(dummy)
577 c call dpoco(d,maxpoints,npoi1,rcond,help,info)
578 anorm=dlansy('1','U',npoi1,d,maxpoints,help)
579 call dpotrf('U',npoi1,d,maxpoints,info)
580 write (6,*) 'd matrix has been decomposed, time=',
581 c dtime(dummy),dummy(1),dummy(2),info
582 write (6,*) 'start finding condition number'
583 call dpocon('U',npoi1,d,maxpoints,anorm,rcond,help,iwork,info1)
584 write (6,*) 'condition number has been found, time=',
585 c dtime(dummy),dummy(1),dummy(2)
586 write (6,*) rcond,info,info1,anorm
587 write (6,*) 'start finding solution'
588 time1=dtime(dummy)
589 tplarson 1.1 c call dposl(d,maxpoints,npoi1,right)
590 call dpotrs('U',npoi1,1,d,maxpoints,right,maxpoints,info)
591 write (6,*) 'solution has been found, time=',
592 c dtime(dummy),dummy(1),dummy(2)
593 c write (21) (right(i),i=1,npoi1)
594 c close(21)
595 c stop
596 write (6,*) 'rcond,info',rcond,info
597 do 1875,j=1,nrad1
598 do 1870,k=1,nsetp1
599 omega(j,k)=right(ijk(j,k))
600 1870 continue
601 1875 continue
602 do 1877,j=1,nrad1
|
651 tplarson 1.1 c l,n,freq(i),mnlm(i1),itnlm(i1),it1nlm(i1),
652 c s,split(i1),sigma(i1)
653 chisq=chisq+((s-split(i1))*sigma1(i1))**2
654 1895 continue
655 1897 continue
656 write (6,*) 'splittings have been calculated, time=',
657 c dtime(dummy),dummy(1),dummy(2)
658 write (6,*) 'start calculating wiggle measures'
659 time1=dtime(dummy)
660 do 1907,k=1,nsetp1
661 do 1905,j=1,nrad1
662 fro(j,k)=0.0
663 fto(j,k)=0.0
664 c omega(j,k)=sin(tmesh1(k))**2*rmesh1(j)**2
665 c omega(j,k)=rmesh1(j)**3
666 1905 continue
667 1907 continue
668 do 1930,j=1,nrad1
669 do 1920,i=max(2,j-1),min(nrad1-1,j+1)
670 xm1=rmesh1(i-1)-rmesh1(i)
671 xp1=rmesh1(i+1)-rmesh1(i)
672 tplarson 1.1 den=xm1*xp1*(xm1-xp1)
673 if (i.eq.j-1) fji=-xm1
674 if (i.eq.j) fji=xm1-xp1
675 if (i.eq.j+1) fji=xp1
676 do 1910,k=1,nsetp1
677 fro(i,k)=fro(i,k)+omega(j,k)*2*fji/den
678 1910 continue
679 1920 continue
680 1930 continue
681 do 1960,k=1,nsetp1
682 do 1950,i=max(1,k-1),min(nsetp1,k+1)
683 t=tmesh1(i)
684 if (i.gt.1) tm1=tmesh1(i-1)
685 if (i.eq.1) tm1=2*tmesh1(1)-tmesh1(2)
686 if (i.lt.nsetp1) tp1=tmesh1(i+1)
687 if (i.eq.nsetp1) tp1=2*tmesh1(nsetp1)-tmesh1(nsetp1-1)
688 xm1=tm1-t
689 xp1=tp1-t
690 den=xm1*xp1*(xm1-xp1)
691 if (i.eq.k-1) fki=-xm1
692 if (i.eq.k) fki=xm1-xp1
693 tplarson 1.1 if (i.eq.k+1) fki=xp1
694 if (i.eq.1.and.k.eq.2) fki=xp1-xm1
695 if (i.eq.nsetp1.and.k.eq.(nsetp1-1)) fki=xp1-xm1
696 do 1940,j=1,nrad1
697 fto(j,i)=fto(j,i)+omega(j,k)*2*fki/den
698 1940 continue
699 1950 continue
700 1960 continue
701 c calculate wigr and wigt
702 wigr=0.0
703 wigt=0.0
704 do 1980,k=1,nsetp1
705 wt=tweights(k)
706 do 1970,j=1,nrad1
707 wigr=wigr+wt*rweights(j)*fr(j)*fro(j,k)**2
708 wigt=wigt+wt*rweights(j)*ft(j)*fto(j,k)**2
709 1970 continue
710 1980 continue
711 write (6,*) 'wiggle measures have been computed, time=',
712 c dtime(dummy),dummy(1),dummy(2)
713 write (6,*) 'chisq=',chisq
714 tplarson 1.1 write (6,*) 'rms=',sqrt(chisq/nnlm)
715 write (6,*) 'wigr, wigt=',wigr,wigt
716 write (6,*) '%',mur,mut,chisq,sqrt(chisq/nnlm),wigr,wigt,
717 c info,rcond
718 c end do
719 c end do
720 c stop
721 if ((icovar.eq.1).or.(iavkern.ge.1)) then
722 write (6,*) 'start inverting d matrix'
723 time1=dtime(dummy)
724 c call dpodi(d,maxpoints,npoi1,det,1)
725 call dpotri('U',npoi1,d,maxpoints,info)
726 c open (21,file=tmpdir//'saved1',form='unformatted',
727 c c status='unknown',access='direct',recl=8*npoi1)
728 do 2095,i=1,npoi1
729 c do 2092,j=1,npoi1
730 c di(j)=0.0D0
731 c2092 continue
732 c di(i)=1.0D0
733 c call dposl(d,maxpoints,npoi1,di)
734 c write (21,rec=i) (di(j),j=1,npoi1)
735 tplarson 1.1 c write (21,rec=i) (d(j,i),j=1,i),(d(i,j),j=i+1,npoi1)
736 do 2092,j=i+1,npoi1
737 d(j,i)=d(i,j)
738 2092 continue
739 2095 continue
740 c close(21)
741 write (6,*) 'd matrix has been inverted, time=',
742 c dtime(dummy),dummy(1),dummy(2)
743 end if
744 if (icovar.ge.1) then
745 write (6,*) 'start finding variances'
746 c the following loop could be blocked
747 c do 2150,i=1,npoi1
748 c call dgemv('n',npoi1,npoi1,1.0D0,b,maxpoints,d(1,i),1,
749 c c 0.0D0,dib,1)
750 c var(i)=ddot(npoi1,d(1,i),1,dib,1)
751 c2150 continue
752 do 2150,i=1,npoi1,nblock1
753 ix=min(npoi1,i+nblock1-1)
754 nx=ix+1-i
755 call dgemm('n','n',npoi1,nx,npoi1,1.0D0,b,maxpoints,
756 tplarson 1.1 c d(1,i),maxpoints,0.0D0,dibx,maxpoints)
757 do 2140,j=i,ix
758 var(j)=ddot1(npoi1,d(1,j),dibx(1,j-i+1))
759 2140 continue
760 2150 continue
761 do 2190,j=1,nrad1
762 write (26,'(100f11.5)') (sqrt(var(ijk(j,k))),k=1,nsetp1)
763 2190 continue
764 c close(21)
765 write (6,*) 'variances have been found, time=',
766 c dtime(dummy),dummy(1),dummy(2)
767 c else
768 c write (6,*) 'start finding covariance matrix'
769 c time1=dtime(dummy)
770 c do 2230,i=1,npoi1
771 c do 2220,j=1,npoi1
772 c d(j,i)=ddot1(npoi1,b(1,j),d1(1,i))
773 c2220 continue
774 c2230 continue
775 c do 2250,i=1,npoi1
776 c do 2240,j=i,npoi1
777 tplarson 1.1 c cov(j,i)=ddot1(npoi1,d1(1,j),d(1,i))
778 c2240 continue
779 c2250 continue
780 c do 2270,i=1,npoi1
781 c do 2260,j=1,i-1
782 c cov(j,i)=cov(i,j)
783 c2260 continue
784 c2270 continue
785 c do 2280,i=1,nrad1
786 c write (6,*) rmesh1(i)
787 c2280 continue
788 c do 2290,j=1,nrad1
789 c write (6,'(100f9.5)') (sqrt(cov(ijk(j,k),ijk(j,k))),
790 c c k=1,nsetp1)
791 c2290 continue
792 c write (6,*) 'covariance matrix has been found, time=',
793 c c dtime(dummy),dummy(1),dummy(2)
794 end if
795 if (iavkern.eq.0) then
796 write (6,*) 'skipping calculation of averaging kernels'
797 else
798 tplarson 1.1 write (6,*) 'start finding averaging kernels'
799 time1=dtime(dummy)
800 write (6,*) nt2,nr2
801 c open (21,file=tmpdir//'saved1',form='unformatted',
802 c c status='old',access='direct',recl=8*npoi1)
803 c moved up
804 c open (27,file='avkern.2d',form='unformatted')
805 c open (28,file='coeff',form='unformatted')
806 do 2300,j1=1,nr2
807 x2(j1)=radmesh(nsr2*(j1-1)+1)
808 2300 continue
809 pi=4.0D0*atan(1.0D0)
810 pi2=pi/2.0D0
811 do 2301,k1=1,nt2
812 x1(k1)=pi2*(k1-1)/(nt2-1.)
813 2301 end do
814 write (6,*) nrad1,nsr,nsetp1,nst
815 if (nsr.ne.0) then
816 i=0
817 do 2303,j=1,nrad1,nsr
818 do 2302,k=1,nsetp1,nst
819 tplarson 1.1 i=i+1
820 jav(i)=j
821 kav(i)=k
822 2302 continue
823 2303 continue
824 nav=i
825 end if
826 c do 2380,j=1,nrad1,nsr
827 c do 2370,k=1,nsetp1,nst
828 do 2380,iav=1,nav
829 j=jav(iav)
830 k=kav(iav)
831 ipoi=ijk(j,k)
832 c read (21,rec=ipoi) (di(j1),j1=1,npoi1)
833 do 2305,is=1,nnlm
834 coeff(is)=0.0D0
835 2305 continue
836 do 2330,i1=1,npoints
837 i1a=max(1,i1-nident)
838 k1=mod(i1-1,nsetp1)+1
839 j1=(i1-1)/nsetp1+1
840 tplarson 1.1 do 2310,i=1,nmodes
841 f1x(i)=f1(i,j1)
842 f2x(i)=f2(i,j1)
843 2310 continue
844 do 2315,i=1,nlm
845 g1x(i)=g1(i,k1)
846 g2x(i)=g2(i,k1)
847 2315 continue
848 c d1i=d1(i1a,ipoi)
849 c d1i=di(i1a)
850 d1i=d(i1a,ipoi)
851 do 2320,is=1,nnlm
852 i=iln(is)
853 ilmi=ilm(is)
854 coeff(is)=coeff(is)+d1i*sigma1(is)*
855 c (f1x(i)*g1x(ilmi)+f2x(i)*g2x(ilmi))
856 2320 continue
857 2330 continue
858 c s1=0.0D0
859 c s2=0.0D0
860 c s3=0.0D0
861 tplarson 1.1 vari=0.0
862 do 2350,is=1,nnlm
863 coeff(is)=coeff(is)*sigma1(is)
864 c s1=s1+coeff(is)*split(is)
865 c s2=s2+mnlm(is)*coeff(is)
866 c s3=s3+coeff(is)
867 vari=vari+(coeff(is)*sigma(is))**2
868 2350 continue
869 c write (6,*) '#',j,k
870 c write (6,*) s1,right(ipoi),s1-right(ipoi),s2,s3
871 c call set2davk(nnlm,isplit,msplit,coeff,avkern,nsr2,nt2)
872 call set2davk(coeff,avkern,nsr2,nt2)
873 if (iavkern.eq.2) then
874 write (27) ((avkern(k1,j1),k1=1,nt2),j1=1,nr2)
875 end if
876 c call fitmn(x1,x2,nt2,nr2,avkern,yy2,mjta,4,
877 c c real(tmesh1(k)),real(rmesh1(j)))
878 write (24,*) mur,mut,real(rmesh1(j)),real(tmesh1(k)),
879 c right(ipoi),sqrt(vari),
880 c mjta(3),mjta(1),mjta(4),mjta(2)
881 if (icoeff.eq.1) then
882 tplarson 1.1 write (28) (real(coeff(is)),is=1,nnlm)
883 end if
884 c end if
885 c2370 continue
886 2380 continue
887 c close(27)
888 c close(28)
889 c close(21)
890 write (6,*) 'averaging kernels set up, time=',
891 c dtime(dummy),dummy(1),dummy(2)
892 end if
893 5000 continue
894 if (iavkern.ne.0) then
895 close (24)
896 end if
897 close (25)
898 if (icovar.eq.26) then
899 close (26)
900 end if
901 close (27)
902 if (icoeff.eq.28) then
903 tplarson 1.1 close (28)
904 end if
905 open (20,file='status.log',form='formatted')
906 write (20,*) 0
907 close (20)
908
909 write (6,*) 'total time used:',etime(dummy)
910 write (6,*) 'user time: ',dummy(1)
911 write (6,*) 'system time: ',dummy(2)
912 end
913
|
Return to
inv2d.f
CVS log
Up to [Development] / JSOC / proj / globalhs / apps