JSOC/proj/globalhs/apps/inv2d.f - annotate - 1.3

Return to inv2d.f CVS log

Up to [Development] / JSOC / proj / globalhs / apps

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=maxr1maxt1,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 character80 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(3maxpoints) 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 tplarson 1.3 70 open (20,file='cvstag.log',form='formatted') 71 write (20,*) CVSTAG 72 close (20) 73
74 tplarson 1.1 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 tplarson 1.1 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.0D0atan(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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 c dtime(dummy),dummy(1),dummy(2) 285 end if 286 call dsetint(1,nsetp1,tmesh1,tweights) 287 npoints=nrad1nsetp1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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)/(nnlmnpoints),(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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 c just a try 453 c ft(i)=1.0 454 1400 continue 455 do 1405,i=1,nrad1 456 fwr(i)=murrweights(i)fr(i) 457 c fwr(i)=murrweights(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=xm1xp1(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+fjifkifwr(i)/den2 473 tplarson 1.1 1410 continue 474 c factor of 4=2^2 since second derivative is 2a 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)=muttweights(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=2tmesh1(1)-tmesh1(2) 489 if (i.lt.nsetp1) tp1=tmesh1(i+1) 490 if (i.eq.nsetp1) tp1=2tmesh1(nsetp1)-tmesh1(nsetp1-1) 491 xm1=tm1-t 492 xp1=tp1-t 493 den=xm1xp1(xm1-xp1) 494 tplarson 1.1 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+fjifkifwt(i)/den2 505 c write (6,) i,j,k,ckj 506 1510 continue 507 c factor of 4=2^2 since second derivative is 2a 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.2 write (25,'(101f12.4)') (right(ijk(j,k)),k=1,nsetp1)
609 tplarson 1.1 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 tplarson 1.2 open(10,file='splittings.out',form='formatted')
629 tplarson 1.1 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 tplarson 1.1 s=s+g1(index,k)f1o(k)+g2(index,k)f2o(k) 651 1893 continue 652 l=lmode(i) 653 n=nmode(i)
654 tplarson 1.2 write (10,'(i4,i4,f8.2,i4,2i3,3f12.4)') 655 c write (6,'(i4,i4,f8.2,i4,2i3,3f12.4)')
656 tplarson 1.1 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))2rmesh1(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 tplarson 1.1 den=xm1xp1(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)2fji/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=2tmesh1(1)-tmesh1(2) 691 if (i.lt.nsetp1) tp1=tmesh1(i+1) 692 if (i.eq.nsetp1) tp1=2tmesh1(nsetp1)-tmesh1(nsetp1-1) 693 xm1=tm1-t 694 xp1=tp1-t 695 den=xm1xp1(xm1-xp1) 696 if (i.eq.k-1) fki=-xm1 697 if (i.eq.k) fki=xm1-xp1 698 tplarson 1.1 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)2fki/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+wtrweights(j)fr(j)fro(j,k)2 713 wigt=wigt+wtrweights(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 tplarson 1.1 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=8npoi1) 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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=8npoi1) 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.0D0atan(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 tplarson 1.1 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 tplarson 1.1 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)+d1isigma1(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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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

Powered by
ViewCVS 0.9.4