 |
Millepede-II
V04-00-00
|
|
Millepede-II
V04-00-00
|
00001 00002 ! Code converted using TO_F90 by Alan Miller 00003 ! Date: 2012-03-16 Time: 11:06:00 00004 00087 00194 00348 00350 PROGRAM mptwo 00351 USE mpmod 00352 USE mpdalc 00353 USE mptest1, ONLY: nplan,del,dvd 00354 USE mptest2, ONLY: nlyr,nmx,nmy,sdevx,sdevy 00355 00356 IMPLICIT NONE 00357 REAL :: andf 00358 REAL :: c2ndf 00359 REAL :: deltat 00360 REAL :: diff 00361 REAL :: err 00362 REAL :: gbu 00363 REAL :: gmati 00364 REAL :: rej 00365 REAL :: rloop1 00366 REAL :: rloop2 00367 REAL :: rstext 00368 REAL :: secnd 00369 REAL :: rst 00370 REAL ::rstp 00371 REAL, DIMENSION(2) :: ta 00372 INTEGER :: i 00373 INTEGER :: ii 00374 INTEGER :: ix 00375 INTEGER :: ixv 00376 INTEGER :: iy 00377 INTEGER :: k 00378 INTEGER :: kfl 00379 INTEGER :: lun 00380 INTEGER :: minut 00381 INTEGER :: nhour 00382 INTEGER :: nmxy 00383 INTEGER :: nrc 00384 INTEGER :: nsecnd 00385 INTEGER :: ntot 00386 INTEGER :: ntsec 00387 00388 CHARACTER (LEN=24) :: chdate 00389 CHARACTER (LEN=24) :: chost 00390 00391 INTEGER(KIND=large) :: rows 00392 INTEGER(KIND=large) :: cols 00393 00394 DOUBLE PRECISION :: sums(9) 00395 !$ INTEGER :: OMP_GET_NUM_PROCS,OMP_GET_MAX_THREADS 00396 !$ INTEGER :: MXTHRD 00397 !$ INTEGER :: NPROC 00398 00399 SAVE 00400 ! ... 00401 CALL etime(ta,rstp) 00402 CALL fdate(chdate) 00403 00404 ! millepede.log file 00405 lunlog=8 00406 lvllog=1 00407 CALL mvopen(lunlog,'millepede.log') 00408 CALL getenv('HOSTNAME',chost) 00409 IF (chost(1:1) == ' ') CALL getenv('HOST',chost) 00410 WRITE(*,*) '($Rev: 92 $)' 00411 !$ WRITE(*,*) 'using OpenMP (TM)' 00412 !#ifdef __GFORTRAN__ 00413 ! WRITE(*,111) __GNUC__ , __GNUC_MINOR__ , __GNUC_PATCHLEVEL__ 00414 !111 FORMAT(' compiled with gcc ',i0,'.',i0,'.',i0) 00415 !#endif 00416 WRITE(*,*) ' ' 00417 WRITE(*,*) ' < Millepede II-P starting ... ',chdate 00418 WRITE(*,*) ' ',chost 00419 WRITE(*,*) ' ' 00420 00421 WRITE(8,*) ' ' 00422 WRITE(8,*) 'Log-file Millepede II-P ', chdate 00423 WRITE(8,*) ' ', chost 00424 ! CALL MEGDEB ! debug switched on 00425 00426 ! read command line and text files 00427 00428 CALL filetc ! command line and steering file analysis 00429 CALL filetx ! read text files 00430 lvllog=mprint ! export print level 00431 IF (memdbg > 0) printflagalloc=1 ! debug memory management 00432 !$ WRITE(*,*) 00433 !$ NPROC=1 00434 !$ MXTHRD=1 00435 !$ NPROC=OMP_GET_NUM_PROCS() ! number of processors available 00436 !$ CALL OMP_SET_NUM_THREADS(MTHRD) ! set max number of threads to MTHRD 00437 !$ MXTHRD=OMP_GET_MAX_THREADS() ! get max number of threads back 00438 !$ WRITE(*,*) 'Number of processors available: ', NPROC 00439 !$ WRITE(*,*) 'Maximum number of OpenMP threads: ', MXTHRD 00440 !$ WRITE(*,*) 'Number of threads for processing: ', MTHRD 00441 !$ IF (MXREC.GT.0) MTHRDR=1 ! to get allways the same MXREC records 00442 !$ WRITE(*,*) 'Number of threads for reading: ', MTHRDR 00443 !$POMP INST INIT ! start profiling with ompP 00444 IF (ncache < 0) THEN 00445 ncache=25000000*mthrd ! default cache size (100 MB per thread) 00446 ENDIF 00447 rows=6; cols=mthrdr 00448 CALL mpalloc(readBufferInfo,rows,cols,'read buffer header') 00449 ! histogram file 00450 lun=7 00451 CALL mvopen(lun,'millepede.his') 00452 CALL hmplun(lun) ! unit for histograms 00453 CALL gmplun(lun) ! unit for xy data 00454 00455 ! debugging 00456 IF(nrecpr /= 0.OR.nrecp2 /= 0) THEN 00457 CALL mvopen(1,'mpdebug.txt') 00458 END IF 00459 00460 CALL etime(ta,rstext) 00461 times(0)=rstext-rstp ! time for text processing 00462 00463 ! preparation of data sub-arrays 00464 00465 CALL loop1 00466 CALL etime(ta,rloop1) 00467 times(1)=rloop1-rstext ! time for LOOP1 00468 00469 CALL loop2 00470 IF(chicut /= 0.0) THEN 00471 WRITE(8,*) 'Chi square cut equiv 3 st.dev applied ...' 00472 WRITE(8,*) ' in first iteration with factor',chicut 00473 WRITE(8,*) ' in second iteration with factor',chirem 00474 WRITE(8,*) ' (reduced by sqrt in next iterations)' 00475 END IF 00476 IF(lhuber /= 0) THEN 00477 WRITE(8,*) 'Down-weighting of outliers in', lhuber,' iterations' 00478 WRITE(8,*) 'Cut on downweight fraction',dwcut 00479 END IF 00480 00481 CALL etime(ta,rloop2) 00482 times(2)=rloop2-rloop1 ! time for LOOP2 00483 00484 ! use different solution methods 00485 00486 CALL mstart('Iteration') ! Solution module starting 00487 00488 ! CALL XEXECS ! all methods 00489 00490 CALL xloopn ! all methods 00491 00492 ! IF(NREJEC(1)+NREJEC(2)+NREJEC(3).NE.0) THEN 00493 ! WRITE(*,*) 'Data rejected in last loop: ',NREJEC(1), 00494 ! + ' (Ndf=0) ',NREJEC(2),' (huge) ',NREJEC(3),' (large)' 00495 ! END IF 00496 00497 00498 00499 ! ------------------------------------------------------------------ 00500 00501 IF(nloopn > 2.AND.nhistp /= 0) THEN ! last iteration 00502 CALL hmprnt(3) ! scaled residual of single measurement (with global deriv.) 00503 CALL hmprnt(12) ! scaled residual of single measurement (no global deriv.) 00504 CALL hmprnt(4) ! chi^2/Ndf 00505 END IF 00506 IF(nloopn > 2) THEN 00507 CALL hmpwrt(3) 00508 CALL hmpwrt(12) 00509 CALL hmpwrt(4) 00510 CALL gmpwrt(4) ! location, dispersion (res.) as a function of record nr 00511 IF (nloopn <= lfitnp) THEN 00512 CALL hmpwrt(13) 00513 CALL hmpwrt(14) 00514 CALL gmpwrt(5) 00515 END IF 00516 END IF 00517 IF(nhistp /= 0) THEN 00518 CALL gmprnt(1) 00519 CALL gmprnt(2) 00520 END IF 00521 CALL gmpwrt(1) ! output of xy data 00522 CALL gmpwrt(2) ! output of xy data 00523 ! 'track quality' per binary file 00524 IF (nfilb > 1) THEN 00525 CALL gmpdef(6,1,'log10(#records) vs file number') 00526 CALL gmpdef(7,1,'final rejection fraction vs file number') 00527 CALL gmpdef(8,1, & 00528 'final <Chi^2/Ndf> from accepted local fits vs file number') 00529 CALL gmpdef(9,1, '<Ndf> from accepted local fits vs file number') 00530 00531 DO i=1,nfilb 00532 kfl=kfd(2,i) 00533 nrc=-kfd(1,i) 00534 IF (nrc > 0) THEN 00535 rej=FLOAT(nrc-jfd(kfl))/FLOAT(nrc) 00536 CALL gmpxy(6,FLOAT(kfl),ALOG10(FLOAT(nrc))) ! log10(#records) vs file 00537 CALL gmpxy(7,FLOAT(kfl),rej) ! rejection fraction vs file 00538 END IF 00539 IF (jfd(kfl) > 0) THEN 00540 c2ndf=cfd(kfl)/FLOAT(jfd(kfl)) 00541 CALL gmpxy(8,FLOAT(kfl),c2ndf) ! <Chi2/NDF> vs file 00542 andf=FLOAT(dfd(kfl))/FLOAT(jfd(kfl)) 00543 CALL gmpxy(9,FLOAT(kfl),andf) ! <NDF> vs file 00544 END IF 00545 END DO 00546 IF(nhistp /= 0) THEN 00547 CALL gmprnt(6) 00548 CALL gmprnt(7) 00549 CALL gmprnt(8) 00550 CALL gmprnt(9) 00551 END IF 00552 CALL gmpwrt(6) ! output of xy data 00553 CALL gmpwrt(7) ! output of xy data 00554 CALL gmpwrt(8) ! output of xy data 00555 CALL gmpwrt(9) ! output of xy data 00556 END IF 00557 00558 IF(ictest == 1) THEN 00559 WRITE(*,*) ' ' 00560 WRITE(*,*) 'Misalignment test wire chamber' 00561 WRITE(*,*) ' ' 00562 00563 CALL hmpdef( 9,-0.0015,+0.0015,'True - fitted displacement') 00564 CALL hmpdef(10,-0.0015,+0.0015,'True - fitted Vdrift') 00565 DO i=1,4 00566 sums(i)=0.0D0 00567 END DO 00568 DO i=1,nplan 00569 diff=SNGL(-del(i)-globalParameter(i)) 00570 sums(1)=sums(1)+diff 00571 sums(2)=sums(2)+diff*diff 00572 diff=SNGL(-dvd(i)-globalParameter(100+i)) 00573 sums(3)=sums(3)+diff 00574 sums(4)=sums(4)+diff*diff 00575 END DO 00576 sums(1)=0.01D0*sums(1) 00577 sums(2)=SQRT(0.01D0*sums(2)) 00578 sums(3)=0.01D0*sums(3) 00579 sums(4)=SQRT(0.01D0*sums(4)) 00580 WRITE(*,143) 'Parameters 1 - 100: mean =',sums(1), 'rms =',sums(2) 00581 WRITE(*,143) 'Parameters 101 - 200: mean =',sums(3), 'rms =',sums(4) 00582 143 FORMAT(6X,a28,f9.6,3X,a5,f9.6) 00583 WRITE(*,*) ' ' 00584 WRITE(*,*) ' ' 00585 WRITE(*,*) ' I ' 00586 WRITE(*,*) ' --- ' 00587 DO i=1,100 00588 WRITE(*,102) i,-del(i),globalParameter(i),-del(i)-globalParameter(i), & 00589 -dvd(i),globalParameter(100+i),-dvd(i)-globalParameter(100+i) 00590 diff=SNGL(-del(i)-globalParameter(i)) 00591 CALL hmpent( 9,diff) 00592 diff=SNGL(-dvd(i)-globalParameter(100+i)) 00593 CALL hmpent(10,diff) 00594 END DO 00595 IF(nhistp /= 0) THEN 00596 CALL hmprnt( 9) 00597 CALL hmprnt(10) 00598 END IF 00599 CALL hmpwrt( 9) 00600 CALL hmpwrt(10) 00601 END IF 00602 IF(ictest > 1) THEN 00603 WRITE(*,*) ' ' 00604 WRITE(*,*) 'Misalignment test Si tracker' 00605 WRITE(*,*) ' ' 00606 00607 CALL hmpdef( 9,-0.0025,+0.0025,'True - fitted displacement X') 00608 CALL hmpdef(10,-0.025,+0.025,'True - fitted displacement Y') 00609 DO i=1,9 00610 sums(i)=0.0D0 00611 END DO 00612 nmxy=nmx*nmy 00613 ix=0 00614 iy=ntot 00615 DO i=1,nlyr 00616 DO k=1,nmxy 00617 ix=ix+1 00618 diff=SNGL(-sdevx((i-1)*nmxy+k)-globalParameter(ix)) 00619 sums(1)=sums(1)+1.0D0 00620 sums(2)=sums(2)+diff 00621 sums(3)=sums(3)+diff*diff 00622 ixv=globalParLabelIndex(2,ix) 00623 IF (ixv > 0.AND.metsol == 1.OR.metsol == 2) THEN 00624 ii=(ixv*ixv+ixv)/2 00625 gmati=SNGL(globalMatD(ii)) 00626 ERR=SQRT(ABS(gmati)) 00627 diff=diff/ERR 00628 sums(7)=sums(7)+1.0D0 00629 sums(8)=sums(8)+diff 00630 sums(9)=sums(9)+diff*diff 00631 END IF 00632 END DO 00633 IF (MOD(i,3) == 1) THEN 00634 DO k=1,nmxy 00635 iy=iy+1 00636 diff=-SNGL(sdevy((i-1)*nmxy+k)-globalParameter(iy)) 00637 sums(4)=sums(4)+1.0D0 00638 sums(5)=sums(5)+diff 00639 sums(6)=sums(6)+diff*diff 00640 ixv=globalParLabelIndex(2,iy) 00641 IF (ixv > 0.AND.metsol == 1.OR.metsol == 2) THEN 00642 ii=(ixv*ixv+ixv)/2 00643 gmati=SNGL(globalMatD(ii)) 00644 ERR=SQRT(ABS(gmati)) 00645 diff=diff/ERR 00646 sums(7)=sums(7)+1.0D0 00647 sums(8)=sums(8)+diff 00648 sums(9)=sums(9)+diff*diff 00649 END IF 00650 END DO 00651 END IF 00652 END DO 00653 sums(2)=sums(2)/sums(1) 00654 sums(3)=SQRT(sums(3)/sums(1)) 00655 sums(5)=sums(5)/sums(4) 00656 sums(6)=SQRT(sums(6)/sums(4)) 00657 WRITE(*,143) 'Parameters 1 - 500: mean =',sums(2), 'rms =',sums(3) 00658 WRITE(*,143) 'Parameters 501 - 700: mean =',sums(5), 'rms =',sums(6) 00659 IF (sums(7) > 0.5D0) THEN 00660 sums(8)=sums(8)/sums(7) 00661 sums(9)=SQRT(sums(9)/sums(7)) 00662 WRITE(*,143) 'Parameter pulls, all: mean =',sums(8), 'rms =',sums(9) 00663 END IF 00664 WRITE(*,*) ' ' 00665 WRITE(*,*) ' ' 00666 WRITE(*,*) ' I ' 00667 WRITE(*,*) ' --- ' 00668 ix=0 00669 iy=ntot 00670 DO i=1,nlyr 00671 DO k=1,nmxy 00672 ix=ix+1 00673 diff=SNGL(-sdevx((i-1)*nmxy+k)-globalParameter(ix)) 00674 CALL hmpent( 9,diff) 00675 WRITE(*,102) ix,-sdevx((i-1)*nmxy+k),globalParameter(ix),-diff 00676 END DO 00677 END DO 00678 DO i=1,nlyr 00679 IF (MOD(i,3) == 1) THEN 00680 DO k=1,nmxy 00681 iy=iy+1 00682 diff=SNGL(-sdevy((i-1)*nmxy+k)-globalParameter(iy)) 00683 CALL hmpent(10,diff) 00684 WRITE(*,102) iy,-sdevy((i-1)*nmxy+k),globalParameter(iy),-diff 00685 END DO 00686 END IF 00687 END DO 00688 IF(nhistp /= 0) THEN 00689 CALL hmprnt( 9) 00690 CALL hmprnt(10) 00691 END IF 00692 CALL hmpwrt( 9) 00693 CALL hmpwrt(10) 00694 END IF 00695 00696 IF(nrec1+nrec2 > 0) THEN 00697 WRITE(8,*) ' ' 00698 IF(nrec1 > 0) THEN 00699 WRITE(8,*) 'Record',nrec1,' has largest residual:',value1 00700 END IF 00701 IF(nrec2 > 0) THEN 00702 WRITE(8,*) 'Record',nrec2,' has largest Chi^2/Ndf:',value2 00703 END IF 00704 END IF 00705 IF(nrec3 < maxi4) THEN 00706 WRITE(8,*) 'Record',nrec3, ' is first with error (rank deficit/NaN)' 00707 END IF 00708 WRITE(8,*) ' ' 00709 WRITE(8,*) 'In total 2 +',nloopn,' loops through the data files' 00710 IF (mnrsit > 0) THEN 00711 WRITE(8,*) ' ' 00712 WRITE(8,*) 'In total ',mnrsit,' internal MINRES iterations' 00713 END IF 00714 00715 WRITE(8,103) times(0),times(1),times(2),times(4),times(7), & 00716 times(5),times(8),times(3),times(6) 00717 CALL etime(ta,rst) 00718 deltat=rst-rstp 00719 ntsec=nint(deltat) 00720 CALL sechms(deltat,nhour,minut,secnd) 00721 nsecnd=nint(secnd) ! round 00722 WRITE(8,*) 'Total time =',ntsec,' seconds =',nhour,' h',minut, & 00723 ' m',nsecnd,' seconds' 00724 CALL fdate(chdate) 00725 WRITE(8,*) 'end ', chdate 00726 gbu=4.02E-9*FLOAT(maxwordsalloc) ! GB used 00727 WRITE(8,*) ' ' 00728 WRITE(8,105) gbu 00729 00730 ! Rejects ---------------------------------------------------------- 00731 00732 IF(nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) /= 0) THEN 00733 WRITE(8,*) ' ' 00734 WRITE(8,*) 'Data rejected in last iteration: ' 00735 WRITE(8,*) ' ', & 00736 nrejec(0), ' (rank deficit/NaN) ',nrejec(1),' (Ndf=0) ', & 00737 nrejec(2), ' (huge) ',nrejec(3),' (large)' 00738 WRITE(8,*) ' ' 00739 END IF 00740 CALL explfc(8) 00741 00742 WRITE(*,*) ' ' 00743 WRITE(*,*) ' < Millepede II-P ending ... ', chdate ! with exit code',ITEXIT,' >' 00744 WRITE(*,*) ' ' 00745 gbu=4.0E-9*FLOAT(maxwordsalloc) ! GB used 00746 WRITE(*,105) gbu 00747 WRITE(*,*) ' ' 00748 00749 102 FORMAT(2X,i4,2X,3F10.5,2X,3F10.5) 00750 103 FORMAT(' Times [in sec] for text processing',f12.3/ & 00751 ' LOOP1',f12.3/ & 00752 ' LOOP2',f12.3/ & 00753 ' func. value ',f12.3,' *',f4.0/ & 00754 ' func. value, global matrix, solution',f12.3,' *',f4.0/ & 00755 ' new solution',f12.3,' *',f4.0/) 00756 105 FORMAT(' Peak dynamic memory allocation: ',f11.6,' GB') 00757 END PROGRAM mptwo ! Mille 00758 00759 00766 00767 SUBROUTINE solglo(ivgbi) 00768 USE mpmod 00769 00770 IMPLICIT NONE 00771 REAL :: dpa 00772 REAL :: err 00773 REAL :: gcor2 00774 REAL :: par 00775 INTEGER :: iph 00776 INTEGER :: istop 00777 INTEGER :: itgbi 00778 INTEGER :: itgbl 00779 INTEGER :: itn 00780 INTEGER :: itnlim 00781 INTEGER :: nout 00782 00783 INTEGER, INTENT(IN) :: ivgbi 00784 00785 DOUBLE PRECISION :: shift 00786 DOUBLE PRECISION :: rtol 00787 DOUBLE PRECISION :: anorm 00788 DOUBLE PRECISION :: acond 00789 DOUBLE PRECISION :: rnorm 00790 DOUBLE PRECISION :: ynorm 00791 DOUBLE PRECISION :: gmati 00792 DOUBLE PRECISION :: diag 00793 INTEGER(KIND=large) :: ijadd 00794 INTEGER(KIND=large) :: jk 00795 INTEGER(KIND=large) :: ii 00796 LOGICAL :: checka 00797 EXTERNAL avprod, mcsolv, mvsolv 00798 SAVE 00799 DATA iph/0/ 00800 ! ... 00801 IF(iph == 0) THEN 00802 iph=1 00803 WRITE(*,101) 00804 END IF 00805 itgbi=globalParVarToTotal(ivgbi) 00806 itgbl=globalParLabelIndex(1,itgbi) 00807 00808 globalVector=0.0D0 ! reset rhs vector IGVEC 00809 globalVector(ivgbi)=1.0D0 00810 00811 ! NOUT =6 00812 nout =0 00813 itnlim=200 00814 shift =0.0D0 00815 rtol = mrestl ! from steering 00816 checka=.FALSE. 00817 00818 IF(mbandw == 0) THEN ! default preconditioner 00819 CALL minres(nagb,globalVector, & 00820 workspaceD(1), workspaceD(1+nagb), workspaceD(1+2*nagb), & 00821 workspaceD(1+3*nagb),workspaceD(1+4*nagb),workspaceD(1+5*nagb), & 00822 globalCorrections,workspaceMinres, avprod, mcsolv, checka ,.TRUE. , shift, & 00823 nout , itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm) 00824 ELSE IF(mbandw > 0) THEN ! band matrix preconditioner 00825 CALL minres(nagb,globalVector, & 00826 workspaceD(1), workspaceD(1+nagb), workspaceD(1+2*nagb), & 00827 workspaceD(1+3*nagb),workspaceD(1+4*nagb),workspaceD(1+5*nagb), & 00828 globalCorrections,workspaceMinres, avprod, mvsolv, checka ,.TRUE. , shift, & 00829 nout, itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm) 00830 ELSE 00831 CALL minres(nagb,globalVector, & 00832 workspaceD(1), workspaceD(1+nagb), workspaceD(1+2*nagb), & 00833 workspaceD(1+3*nagb),workspaceD(1+4*nagb),workspaceD(1+5*nagb), & 00834 globalCorrections,workspaceMinres, avprod, mvsolv, checka ,.FALSE., shift, & 00835 nout, itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm) 00836 END IF 00837 00838 ! subroutine MINRES( n, b, r1, r2, v, w, w1, w2, x, y, 00839 ! $ AVPROD, Msolve, .TRUE. , .FALSE. , SHIFT, 00840 ! $ NOUT , ITNLIM, rtol, 00841 ! $ ISTOP, ITN, Anorm, Acond, rnorm, ynorm ) 00842 00843 par=SNGL(globalParameter(itgbi)) 00844 dpa=par-globalParStart(itgbi) 00845 gmati=globalCorrections(ivgbi) 00846 ERR=SQRT(ABS(SNGL(gmati))) 00847 IF(gmati < 0.0D0) ERR=-ERR 00848 IF(matsto == 1) THEN ! normal matrix ! ??? 00849 ii=ivgbi 00850 jk=(ii*ii+ii)/2 00851 ELSE IF(matsto == 2) THEN ! sparse matrix 00852 jk=ijadd(ivgbi,ivgbi) 00853 END IF 00854 IF (jk > 0) THEN 00855 diag=globalMatD(jk) 00856 ELSE 00857 diag=DBLE(globalMatF(-jk)) 00858 END IF 00859 gcor2=SNGL(1.0D0-1.0D0/(gmati*diag)) ! global correlation (squared) 00860 WRITE(*,102) itgbl,par,globalParPreSigma(itgbi),dpa,ERR,gcor2,itn 00861 101 FORMAT(1X,' label parameter presigma differ', & 00862 ' Error gcor^2 iit'/ 1X,'---------',2X,5('-----------'),2X,'----') 00863 102 FORMAT(i10,2X,4G12.4,f7.4,i6,i4) 00864 END SUBROUTINE solglo 00865 00867 SUBROUTINE addcst 00868 USE mpmod 00869 00870 IMPLICIT NONE 00871 REAL :: climit 00872 REAL :: factr 00873 REAL :: hugeVal 00874 REAL :: sgm 00875 00876 INTEGER :: i 00877 INTEGER :: icgb 00878 INTEGER :: irhs 00879 INTEGER :: itgbi 00880 INTEGER :: ivgb 00881 INTEGER :: l 00882 INTEGER :: label 00883 INTEGER :: nop 00884 INTEGER :: inone 00885 00886 DOUBLE PRECISION :: rhs 00887 DOUBLE PRECISION :: drhs(4) 00888 INTEGER :: idrh (4) 00889 SAVE 00890 ! ... 00891 nop=0 00892 hugeVal=1.0 ! E6 00893 IF(lenConstraints == 0) RETURN ! no constraints 00894 climit=1.0E-5 ! limit for printout 00895 irhs=0 ! number of values in DRHS(.), to be printed 00896 00897 i=1 00898 DO icgb=1,ncgb 00899 rhs=listConstraints(i )%value ! right hand side 00900 sgm=listConstraints(i+1)%value ! sigma parameter 00901 i=i+2 00902 DO 00903 label=listConstraints(i)%label 00904 factr=listConstraints(i)%value 00905 itgbi=inone(label) ! -> ITGBI= index of parameter label 00906 ivgb =globalParLabelIndex(2,itgbi) ! -> variable-parameter index 00907 00908 IF(icalcm == 1.AND.ivgb > 0) THEN 00909 CALL mupdat(nvgb+icgb,ivgb,DBLE(hugeVal*factr)) ! add to matrix 00910 END IF 00911 00912 rhs=rhs-factr*globalParameter(itgbi) ! reduce residuum 00913 i=i+1 00914 IF(i > lenConstraints) EXIT 00915 IF(listConstraints(i)%label == 0) EXIT 00916 END DO 00917 IF(ABS(rhs) > climit) THEN 00918 irhs=irhs+1 00919 idrh(irhs)=icgb 00920 drhs(irhs)=rhs 00921 nop=1 00922 IF(irhs == 4) THEN 00923 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs) 00924 irhs=0 00925 END IF 00926 END IF 00927 IF(rhs == 0.0) rhs=1.0E-9 00928 globalVector(nvgb+icgb)=rhs*hugeVal 00929 END DO 00930 00931 IF(irhs /= 0) THEN 00932 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs) 00933 END IF 00934 IF(nop == 0) RETURN 00935 WRITE(*,102) ' Constraints: only equation values >', climit,' are printed' 00936 101 FORMAT(' ',4(i4,g11.3)) 00937 102 FORMAT(a,g11.2,a) 00938 END SUBROUTINE addcst 00939 00943 00944 SUBROUTINE feasma 00945 USE mpmod 00946 USE mpdalc 00947 00948 IMPLICIT NONE 00949 REAL :: factr 00950 REAL :: sgm 00951 INTEGER :: i 00952 INTEGER :: icgb 00953 INTEGER :: ij 00954 INTEGER :: itgbi 00955 INTEGER :: ivgb 00956 INTEGER :: j 00957 INTEGER :: l 00958 INTEGER :: label 00959 INTEGER :: nrank 00960 INTEGER :: inone 00961 00962 DOUBLE PRECISION:: rhs 00963 INTEGER(KIND=large):: length 00964 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: vecConstraints 00965 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: matConstraintsT 00966 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: auxVectorD 00967 INTEGER, DIMENSION(:), ALLOCATABLE :: auxVectorI 00968 SAVE 00969 ! ... 00970 IF(lenConstraints == 0) RETURN ! no constraints 00971 00972 length=(ncgb*ncgb+ncgb)/2 00973 CALL mpalloc(matConsProduct, length, 'product matrix of constraints') 00974 matConsProduct=0.0D0 00975 length=ncgb 00976 CALL mpalloc(vecConsResiduals, length, 'residuals of constraints') 00977 CALL mpalloc(vecConstraints,length,'rhs vector') ! double precision vector 00978 CALL mpalloc(auxVectorI,length,'auxiliary array (I)') ! int aux 1 00979 CALL mpalloc(auxVectorD,length,'auxiliary array (D)') ! double aux 1 00980 ! constraint matrix A and product A A^T (A is stored as transposed) 00981 length = ncgb*nvgb 00982 CALL mpalloc(matConstraintsT,length,'transposed matrix of constraints') 00983 matConstraintsT=0.0D0 00984 00985 i=1 00986 DO icgb=1,ncgb 00987 rhs=listConstraints(i )%value ! right hand side 00988 sgm=listConstraints(i+1)%value ! sigma parameter 00989 i=i+2 00990 DO 00991 label=listConstraints(i)%label 00992 factr=listConstraints(i)%value 00993 itgbi=inone(label) ! -> ITGBI= index of parameter label 00994 ivgb =globalParLabelIndex(2,itgbi) ! -> variable-parameter index 00995 IF(ivgb > 0) matConstraintsT(icgb+(ivgb-1)*ncgb)=factr ! matrix element 00996 rhs=rhs-factr*globalParameter(itgbi) ! reduce residuum 00997 i=i+1 00998 IF(i > lenConstraints) EXIT 00999 IF(listConstraints(i)%label == 0) EXIT 01000 END DO 01001 vecConstraints(icgb)=rhs ! constraint discrepancy 01002 END DO 01003 01004 DO l=1,nvgb 01005 ij=0 01006 DO i=1,ncgb 01007 DO j=1,i 01008 ij=ij+1 01009 matConsProduct(ij)=matConsProduct(ij)+matConstraintsT((l-1)*ncgb+i)*matConstraintsT((l-1)*ncgb+j) 01010 END DO 01011 END DO 01012 END DO 01013 01014 ! inversion of product matrix of constraints 01015 CALL sqminv(matConsProduct,vecConstraints,ncgb,nrank, auxVectorD, auxVectorI) 01016 01017 nmiss1=ncgb-nrank 01018 01019 WRITE(*,*) ' ' 01020 WRITE(*,*) 'Rank of product matrix of constraints is',nrank, & 01021 ' for',ncgb,' constraint equations' 01022 WRITE(8,*) 'Rank of product matrix of constraints is',nrank, & 01023 ' for',ncgb,' constraint equations' 01024 IF(nrank < ncgb) THEN 01025 WRITE(*,*) 'Warning: insufficient constraint equations!' 01026 WRITE(8,*) 'Warning: insufficient constraint equations!' 01027 IF (iforce == 0) THEN 01028 isubit=1 01029 WRITE(*,*) ' --> enforcing SUBITO mode' 01030 WRITE(8,*) ' --> enforcing SUBITO mode' 01031 END IF 01032 END IF 01033 01034 CALL mpdealloc(matConstraintsT) 01035 CALL mpdealloc(auxVectorD) 01036 CALL mpdealloc(auxVectorI) 01037 CALL mpdealloc(vecConstraints) 01038 01039 RETURN 01040 END SUBROUTINE feasma ! matrix for feasible solution 01041 01049 SUBROUTINE feasib(concut,iact) 01050 USE mpmod 01051 USE mpdalc 01052 01053 IMPLICIT NONE 01054 REAL :: factr 01055 REAL :: sgm 01056 INTEGER :: i 01057 INTEGER :: icgb 01058 INTEGER :: iter 01059 INTEGER :: itgbi 01060 INTEGER :: ivgb 01061 INTEGER :: label 01062 INTEGER :: inone 01063 01064 REAL, INTENT(IN) :: concut 01065 INTEGER, INTENT(OUT) :: iact 01066 01067 DOUBLE PRECISION :: rhs 01068 DOUBLE PRECISION ::sum1 01069 DOUBLE PRECISION ::sum2 01070 DOUBLE PRECISION ::sum3 01071 INTEGER(KIND=large) :: length 01072 01073 ! ! acts on subarray "0" 01074 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: auxVectorD 01075 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: vecCorrections 01076 SAVE 01077 01078 iact=0 01079 IF(lenConstraints == 0) RETURN ! no constraints 01080 01081 DO iter=1,2 01082 vecConsResiduals=0.0D0 01083 01084 ! calculate right constraint equation discrepancies 01085 i=1 01086 DO icgb=1,ncgb 01087 rhs=listConstraints(i )%value ! right hand side 01088 sgm=listConstraints(i+1)%value ! sigma parameter 01089 i=i+2 01090 DO 01091 label=listConstraints(i)%label 01092 factr=listConstraints(i)%value 01093 itgbi=inone(label) ! -> ITGBI= index of parameter label 01094 rhs=rhs-factr*globalParameter(itgbi) ! reduce residuum 01095 i=i+1 01096 IF(i > lenConstraints) EXIT 01097 IF(listConstraints(i)%label == 0) EXIT 01098 ENDDO 01099 vecConsResiduals(icgb)=rhs ! constraint discrepancy 01100 END DO 01101 01102 ! constraint equation discrepancies ------------------------------- 01103 01104 sum1=0.0D0 01105 sum2=0.0D0 01106 sum3=0.0D0 01107 DO icgb=1,ncgb 01108 sum1=sum1+vecConsResiduals(icgb)**2 01109 sum2=sum2+ABS(vecConsResiduals(icgb)) 01110 sum3=MAX(sum3,ABS(vecConsResiduals(icgb))) 01111 END DO 01112 sum1=SQRT(sum1/dfloat(ncgb)) 01113 sum2=sum2/dfloat(ncgb) 01114 01115 IF(iter == 1.AND.sum1 < concut) RETURN ! do nothing if correction small 01116 01117 IF(iter == 1.AND.ncgb <= 12) THEN 01118 WRITE(*,*) ' ' 01119 WRITE(*,*) 'Constraint equation discrepancies:' 01120 WRITE(*,101) (icgb,vecConsResiduals(icgb),icgb=1,ncgb) 01121 101 FORMAT(4X,4(i5,g12.4)) 01122 WRITE(*,103) concut 01123 103 FORMAT(10X,' Cut on rms value is',g8.1) 01124 END IF 01125 01126 IF(iact == 0) THEN 01127 WRITE(*,*) ' ' 01128 WRITE(*,*) 'Improve constraints' 01129 END IF 01130 iact=1 01131 01132 WRITE(*,102) iter,sum1,sum2,sum3 01133 102 FORMAT(i6,' rms',g12.4,' avrg_abs',g12.4,' max_abs',g12.4) 01134 01135 length=ncgb 01136 CALL mpalloc(auxVectorD,length,'auxiliary vector (D)') 01137 length=nvgb 01138 CALL mpalloc(vecCorrections,length,'constraint corrections') 01139 vecCorrections=0.0D0 01140 01141 ! multiply inverse matrix and constraint vector 01142 CALL dbsvx(matConsProduct,vecConsResiduals,auxVectorD,ncgb) 01143 01144 i=1 01145 DO icgb=1,ncgb 01146 rhs=listConstraints(i )%value ! right hand side 01147 sgm=listConstraints(i+1)%value ! sigma parameter 01148 i=i+2 01149 DO 01150 label=listConstraints(i)%label 01151 factr=listConstraints(i)%value 01152 itgbi=inone(label) ! -> ITGBI= index of parameter label 01153 ivgb =globalParLabelIndex(2,itgbi) ! -> variable-parameter index 01154 IF(ivgb > 0) THEN 01155 vecCorrections(ivgb)=vecCorrections(ivgb)+auxVectorD(icgb)*factr 01156 END IF 01157 i=i+1 01158 IF(i > lenConstraints) EXIT 01159 IF(listConstraints(i)%label == 0) EXIT 01160 ENDDO 01161 END DO 01162 01163 DO i=1,nvgb ! add corrections 01164 itgbi=globalParVarToTotal(i) 01165 globalParameter(itgbi)=globalParameter(itgbi)+vecCorrections(i) 01166 END DO 01167 01168 CALL mpdealloc(vecCorrections) 01169 CALL mpdealloc(auxVectorD) 01170 01171 END DO ! iteration 1 and 2 01172 01173 END SUBROUTINE feasib ! make parameters feasible 01174 01207 SUBROUTINE peread(more) 01208 USE mpmod 01209 01210 IMPLICIT NONE 01211 INTEGER :: i 01212 INTEGER :: iact 01213 INTEGER :: ierrc 01214 INTEGER :: ierrf 01215 INTEGER :: inder 01216 INTEGER :: ioffp 01217 INTEGER :: ithr 01218 INTEGER :: jfile 01219 INTEGER :: jrec 01220 INTEGER :: k 01221 INTEGER :: kfile 01222 INTEGER :: l 01223 INTEGER :: lun 01224 INTEGER :: mpri 01225 INTEGER :: n 01226 INTEGER :: nact 01227 INTEGER :: nbuf 01228 INTEGER :: ndata 01229 INTEGER :: noff 01230 INTEGER :: npointer 01231 INTEGER :: npri 01232 INTEGER :: nr 01233 INTEGER :: nrc 01234 INTEGER :: nrd 01235 INTEGER :: nrpr 01236 INTEGER :: nthr 01237 INTEGER :: ntot 01238 INTEGER :: maxRecordSize 01239 INTEGER :: maxRecordFile 01240 01241 INTEGER, INTENT(OUT) :: more 01242 01243 LOGICAL :: lprint 01244 LOGICAL :: floop 01245 LOGICAL :: eof 01246 DOUBLE PRECISION :: ds0 01247 DOUBLE PRECISION :: ds1 01248 DOUBLE PRECISION :: ds2 01249 DOUBLE PRECISION :: dw 01250 !$ INTEGER :: OMP_GET_THREAD_NUM 01251 SAVE 01252 01253 inder(i)=readBufferDataI(i) 01254 01255 DATA lprint/.TRUE./ 01256 DATA floop/.TRUE./ 01257 DATA npri / 0 /, mpri / 1000 / 01258 ! ... 01259 IF(ifile == 0) THEN ! start/restart 01260 nrec=0 01261 ntot=0 01262 sumRecords=0 01263 skippedRecords=0 01264 numBlocks=0 01265 minRecordsInBlock=size(readBufferDataI) 01266 maxRecordsInBlock=0 01267 readBufferInfo=0 ! reset management info 01268 nrpr=1 01269 nthr=mthrdr 01270 nact=0 ! active threads (have something still to read) 01271 DO k=1,nthr 01272 IF (ifile < nfilb) THEN 01273 ifile=ifile+1 01274 readBufferInfo(1,k)=ifile 01275 readBufferInfo(2,k)=nact 01276 nact=nact+1 01277 END IF 01278 END DO 01279 END IF 01280 nPointer=size(readBufferPointer)/nact 01281 nData=size(readBufferDataI)/nact 01282 more=-1 01283 DO k=1,nthr 01284 iact=readBufferInfo(2,k) 01285 readBufferInfo(4,k)=0 ! reset counter 01286 readBufferInfo(5,k)=iact*nData ! reset offset 01287 END DO 01288 numBlocks=numBlocks+1 ! new block 01289 01290 !$OMP PARALLEL & 01291 !$OMP DEFAULT(PRIVATE) & 01292 !$OMP SHARED(readBufferInfo,readBufferPointer,readBufferDataI,readBufferDataF, & 01293 !$OMP nPointer,nData,skippedRecords,ndimbuf,NTHR,NFILF,FLOOP, & 01294 !$OMP IFD,KFD,IFILE,NFILB,WFD,XFD) & 01295 !$OMP NUM_THREADS(NTHR) 01296 01297 ithr=1 01298 !$ ITHR=OMP_GET_THREAD_NUM()+1 ! thread number 01299 jfile=readBufferInfo(1,ithr) ! file index 01300 iact =readBufferInfo(2,ithr) ! active thread number 01301 jrec =readBufferInfo(3,ithr) ! records read 01302 ioffp=iact*nPointer 01303 01304 files: DO WHILE (jfile > 0) 01305 kfile=kfd(2,jfile) 01306 records: DO 01307 nbuf=readBufferInfo(4,ithr)+1 01308 noff=readBufferInfo(5,ithr)+2 ! 2 header words per record 01309 nr=ndimbuf 01310 IF(kfile <= nfilf) THEN ! Fortran file 01311 lun=kfile+10 01312 READ(lun,IOSTAT=ierrf) n,(readBufferDataF(noff+i),i=1,min(n/2,nr)),& 01313 (readBufferDataI(noff+i),i=1,min(n/2,nr)) 01314 nr=n/2 01315 IF (ierrf < 0) REWIND lun ! end-of-file 01316 eof=(ierrf /= 0) 01317 ELSE ! C file 01318 lun=kfile-nfilf 01319 #ifdef READ_C_FILES 01320 CALL readc(readBufferDataF(noff+1),readBufferDataI(noff+1),nr,lun,ierrc) 01321 n=nr+nr 01322 #else 01323 ierrc=0 01324 #endif 01325 eof=(ierrc <= 0.AND.ierrc /= -4) ! allow buffer overruns -> skip record 01326 END IF 01327 IF(eof) EXIT records ! end-of-files or error 01328 01329 jrec=jrec+1 01330 readBufferInfo(3,ithr)=jrec 01331 IF(floop) THEN 01332 xfd(jfile)=max(xfd(jfile),n) 01333 IF(ithr == 1) THEN 01334 CALL hmplnt(1,n) 01335 IF(inder(noff+1) /= 0) CALL hmpent(8,FLOAT(inder(noff+1))) 01336 END IF 01337 END IF 01338 01339 IF (nr <= ndimbuf) THEN 01340 readBufferInfo(4,ithr)=nbuf 01341 readBufferInfo(5,ithr)=noff+nr 01342 01343 readBufferPointer(ioffp+nbuf)=noff ! pointer to start of buffer 01344 readBufferDataI(noff )=noff+nr ! pointer to end of buffer 01345 readBufferDataI(noff-1)=ifd(kfile)+jrec ! global record number (available with LOOP2) 01346 readBufferDataF(noff )=real(kfile) ! file number 01347 readBufferDataF(noff-1)=wfd(kfile) ! weight 01348 01349 IF ((noff+nr+2+ndimbuf >= nData*(iact+1)).OR.(nbuf >= nPointer)) EXIT files ! buffer full 01350 ELSE 01351 !$OMP ATOMIC 01352 skippedRecords=skippedRecords+1 01353 CYCLE records 01354 END IF 01355 01356 END DO records 01357 01358 readBufferInfo(1,ithr)=-jfile ! flag eof 01359 IF (kfd(1,jfile) == 0) THEN 01360 PRINT *, 'PEREAD: file ', kfile, 'read the first time, found',jrec,' records' 01361 kfd(1,jfile)=-jrec 01362 END IF 01363 ! take next file 01364 !$OMP CRITICAL 01365 IF (ifile < nfilb) THEN 01366 ifile=ifile+1 01367 jrec=0 01368 readBufferInfo(1,ithr)=ifile 01369 readBufferInfo(3,ithr)=jrec 01370 END IF 01371 !$OMP END CRITICAL 01372 jfile=readBufferInfo(1,ithr) 01373 01374 END DO files 01375 !$OMP END PARALLEL 01376 ! compress pointers 01377 nrd=readBufferInfo(4,1) ! buffers from 1 .thread 01378 DO k=2,nthr 01379 iact =readBufferInfo(2,k) 01380 ioffp=iact*nPointer 01381 nbuf=readBufferInfo(4,k) 01382 DO l=1,nbuf 01383 readBufferPointer(nrd+l)=readBufferPointer(ioffp+l) 01384 END DO 01385 nrd=nrd+nbuf 01386 END DO 01387 01388 more=0 01389 DO k=1,nthr 01390 jfile=readBufferInfo(1,k) 01391 IF (jfile > 0) THEN ! no eof yet 01392 readBufferInfo(2,k)=more 01393 more=more+1 01394 ELSE 01395 ! no more files, thread retires 01396 readBufferInfo(1,k)=0 01397 readBufferInfo(2,k)=-1 01398 readBufferInfo(3,k)=0 01399 END IF 01400 END DO 01401 ! record limit ? 01402 IF (mxrec > 0.AND.(ntot+nrd) >= mxrec) THEN 01403 nrd=mxrec-ntot 01404 more=-1 01405 DO k=1,nthr 01406 jfile=readBufferInfo(1,k) 01407 IF (jfile > 0) THEN ! rewind files 01408 nrc=readBufferInfo(3,k) 01409 IF (kfd(1,jfile) == 0) kfd(1,jfile)=-nrc 01410 IF (kfile <= nfilf) THEN 01411 lun=kfile+10 01412 REWIND lun 01413 ELSE 01414 lun=kfile-nfilf 01415 #ifdef READ_C_FILES 01416 CALL resetc(lun) 01417 #endif 01418 END IF 01419 END IF 01420 END DO 01421 END IF 01422 01423 ntot=ntot+nrd 01424 nrec=ntot 01425 numReadbuffer=nrd 01426 01427 sumRecords=sumRecords+nrd 01428 minRecordsInBlock=MIN(minRecordsInBlock,nrd) 01429 maxRecordsInBlock=MAX(maxRecordsInBlock,nrd) 01430 01431 DO WHILE (nloopn == 0.AND.ntot >= nrpr) 01432 WRITE(*,*) ' Record ',nrpr 01433 IF (nrpr < 100000) THEN 01434 nrpr=nrpr*10 01435 ELSE 01436 nrpr=nrpr+100000 01437 END IF 01438 END DO 01439 01440 IF (ncache > 0.AND.nloopn <= 1.AND. npri < mpri.AND.mprint > 1) THEN 01441 npri=npri+1 01442 IF (npri == 1) WRITE(*,100) 01443 WRITE(*,101) nrec, nrd, more ,ifile 01444 100 FORMAT(/' PeRead records active file' & 01445 /' total block threads number') 01446 101 FORMAT(' PeRead',4I10) 01447 END IF 01448 01449 IF (more <= 0) THEN 01450 ifile=0 01451 IF (floop) THEN 01452 ! check for file weights 01453 ds0=0.0D0 01454 ds1=0.0D0 01455 ds2=0.0D0 01456 maxRecordSize=0 01457 maxRecordFile=0 01458 DO k=1,nfilb 01459 IF (xfd(k) > maxRecordSize) THEN 01460 maxRecordSize=xfd(k) 01461 maxRecordFile=k 01462 END IF 01463 dw=dfloat(-kfd(1,k)) 01464 IF (wfd(k) /= 1.0) nfilw=nfilw+1 01465 ds0=ds0+dw 01466 ds1=ds1+dw*DBLE(wfd(k)) 01467 ds2=ds2+dw*DBLE(wfd(k)**2) 01468 END DO 01469 PRINT *, 'PEREAD: file ', maxRecordFile, 'with max record size ', maxRecordSize 01470 IF (nfilw > 0.AND.ds0 > 0.0D0) THEN 01471 ds1=ds1/ds0 01472 ds2=ds2/ds0-ds1*ds1 01473 DO lun=6,lunlog,2 01474 WRITE(lun,177) nfilw,SNGL(ds1),SNGL(ds2) 01475 177 FORMAT(/' !!!!!',i4,' weighted binary files', & 01476 /' !!!!! mean, variance of weights =',2G12.4) 01477 END DO 01478 END IF 01479 ! integrate record numbers 01480 DO k=2,nfilb 01481 ifd(k)=ifd(k-1)-kfd(1,k-1) 01482 END DO 01483 ! sort 01484 IF (nthr > 1) CALL sort2k(kfd,nfilb) 01485 IF (skippedRecords > 0) THEN 01486 PRINT *, 'PEREAD skipped records: ', skippedRecords 01487 ndimbuf=maxRecordSize/2 ! adjust buffer size 01488 END IF 01489 END IF 01490 lprint=.FALSE. 01491 floop=.FALSE. 01492 IF (ncache > 0.AND.nloopn <= 1.AND.mprint > 0) & 01493 WRITE(*,179) numBlocks, sumRecords, minRecordsInBlock, maxRecordsInBlock 01494 179 FORMAT(/' Read cache usage (#blocks, #records, ', & 01495 'min,max records/block'/17X,4I10) 01496 END IF 01497 RETURN 01498 01499 END SUBROUTINE peread 01500 01508 SUBROUTINE peprep(mode) 01509 USE mpmod 01510 01511 IMPLICIT NONE 01512 01513 INTEGER, INTENT(IN) :: mode 01514 01515 INTEGER :: ibuf 01516 INTEGER :: ichunk 01517 INTEGER :: iproc 01518 INTEGER :: isfrst 01519 INTEGER :: islast 01520 INTEGER :: ist 01521 INTEGER :: j 01522 INTEGER :: ja 01523 INTEGER :: jb 01524 INTEGER :: jsp 01525 INTEGER :: nst 01526 INTEGER :: inone 01527 !$ INTEGER :: OMP_GET_THREAD_NUM 01528 01529 01530 isfrst(ibuf)=readBufferPointer(ibuf)+1 01531 islast(ibuf)=readBufferDataI(readBufferPointer(ibuf)) 01532 01533 IF (mode > 0) THEN 01534 ichunk=MIN((numReadBuffer+mthrd-1)/mthrd/32+1,256) 01535 ! parallelize record loop 01536 !$OMP PARALLEL DO & 01537 !$OMP DEFAULT(PRIVATE) & 01538 !$OMP SHARED(numReadBuffer,readBufferPointer,readBufferDataI,readBufferDataF,ICHUNK) & 01539 !$OMP SCHEDULE(DYNAMIC,ICHUNK) 01540 DO ibuf=1,numReadBuffer ! buffer for current record 01541 iproc=0 01542 !$ IPROC=OMP_GET_THREAD_NUM() ! thread number 01543 ist=isfrst(ibuf) 01544 nst=islast(ibuf) 01545 DO ! loop over measurements 01546 CALL isjajb(nst,ist,ja,jb,jsp) 01547 IF(jb == 0) EXIT 01548 DO j=1,ist-jb 01549 readBufferDataI(jb+j)=inone( readBufferDataI(jb+j) ) ! translate to index 01550 END DO 01551 END DO 01552 END DO 01553 !$OMP END PARALLEL DO 01554 END IF 01555 01556 !$POMP INST BEGIN(peprep) 01557 IF (mode <= 0) THEN 01558 DO ibuf=1,numReadBuffer ! buffer for current record 01559 ist=isfrst(ibuf) 01560 nst=islast(ibuf) 01561 DO ! loop over measurements 01562 CALL isjajb(nst,ist,ja,jb,jsp) 01563 IF(jb == 0) EXIT 01564 DO j=1,ist-jb 01565 readBufferDataI(jb+j)=inone( readBufferDataI(jb+j) ) ! translate to index 01566 END DO 01567 END DO 01568 END DO 01569 END IF 01570 !$POMP INST END(peprep) 01571 01572 END SUBROUTINE peprep 01573 01596 SUBROUTINE isjajb(nst,is,ja,jb,jsp) 01597 USE mpmod 01598 01599 IMPLICIT NONE 01600 REAL :: glder 01601 INTEGER :: i 01602 INTEGER :: inder 01603 01604 INTEGER, INTENT(IN) :: nst 01605 INTEGER, INTENT(IN OUT) :: is 01606 INTEGER, INTENT(OUT) :: ja 01607 INTEGER, INTENT(OUT) :: jb 01608 INTEGER, INTENT(OUT) :: jsp 01609 SAVE 01610 ! ... 01611 inder(i)=readBufferDataI(i) 01612 glder(i)=readBufferDataF(i) 01613 01614 jsp=0 01615 DO 01616 ja=0 01617 jb=0 01618 IF(is >= nst) RETURN 01619 DO 01620 is=is+1 01621 IF(inder(is) == 0) EXIT 01622 END DO 01623 ja=is 01624 DO 01625 is=is+1 01626 IF(inder(is) == 0) EXIT 01627 END DO 01628 jb=is 01629 DO WHILE(inder(is+1) /= 0.AND.is < nst) 01630 is=is+1 01631 END DO 01632 IF(ja+1 /= jb.OR.glder(jb) >= 0.0) EXIT 01633 ! special data 01634 jsp=jb ! pointer to special data 01635 is=is+IFIX(-glder(jb)+0.5) ! skip NSP words 01636 END DO 01637 01638 END SUBROUTINE isjajb 01639 01640 01641 !*********************************************************************** 01642 ! LOOPN ... 01648 01649 SUBROUTINE loopn 01650 USE mpmod 01651 01652 IMPLICIT NONE 01653 REAL :: dsum 01654 REAL :: elmt 01655 REAL :: factrj 01656 REAL :: factrk 01657 REAL :: glder 01658 REAL :: peakd 01659 REAL :: peaki 01660 REAL :: ratae 01661 REAL :: rhs 01662 REAL :: rloop 01663 REAL :: sgm 01664 REAL :: used 01665 REAL :: usei 01666 REAL :: weight 01667 INTEGER :: i 01668 INTEGER :: ia 01669 INTEGER :: ib 01670 INTEGER :: ibuf 01671 INTEGER :: inder 01672 INTEGER :: ioffb 01673 INTEGER :: ipr 01674 INTEGER :: isfrst 01675 INTEGER :: islast 01676 INTEGER :: itgbi 01677 INTEGER :: itgbij 01678 INTEGER :: itgbik 01679 INTEGER :: ivgb 01680 INTEGER :: ivgbij 01681 INTEGER :: ivgbik 01682 INTEGER :: j 01683 INTEGER :: k 01684 INTEGER :: lastit 01685 INTEGER :: lun 01686 INTEGER :: ncrit 01687 INTEGER :: ndfs 01688 INTEGER :: ngras 01689 INTEGER :: nparl 01690 INTEGER :: nr 01691 INTEGER :: nrej 01692 INTEGER:: inone 01693 01694 DOUBLE PRECISION:: adder 01695 DOUBLE PRECISION::funref 01696 DOUBLE PRECISION::dchi2s 01697 DOUBLE PRECISION::sndf 01698 INTEGER(KIND=large):: ii 01699 SAVE 01700 ! ... 01701 isfrst(ibuf)=readBufferPointer(ibuf)+1 01702 islast(ibuf)=readBufferDataI(readBufferPointer(ibuf)) 01703 inder(i)=readBufferDataI(i) 01704 glder(i)=readBufferDataF(i) 01705 ! ----- book and reset --------------------------------------------- 01706 IF(nloopn == 0) THEN ! first call 01707 lastit=-1 01708 iitera=0 01709 END IF 01710 01711 nloopn=nloopn+1 ! increase loop counter 01712 ndfsum=0 01713 sumndf=0.0D0 01714 funref=0.0D0 01715 01716 IF(nloopn == 1) THEN ! book histograms for 1. iteration 01717 CALL gmpdef(1,4,'Function value in iterations') 01718 IF (metsol == 3) THEN ! extend to GMRES, i.e. 4? 01719 CALL gmpdef(2,3,'Number of MINRES steps vs iteration nr') 01720 END IF 01721 CALL hmpdef( 5,0.0,0.0,'Number of degrees of freedom') 01722 CALL hmpdef(11,0.0,0.0,'Number of local parameters') 01723 CALL hmpdef(23,0.0,0.0, 'SQRT of diagonal elements without presigma') 01724 CALL hmpdef(24,0.0,0.0, 'Log10 of off-diagonal elements') 01725 CALL hmpdef(25,0.0,0.0, 'Relative individual pre-sigma') 01726 CALL hmpdef(26,0.0,0.0, 'Relative global pre-sigma') 01727 END IF 01728 01729 01730 CALL hmpdef(3,-prange,prange, & ! book 01731 'Normalized residuals of single (global) measurement') 01732 CALL hmpdef(12,-prange,prange, & ! book 01733 'Normalized residuals of single (local) measurement') 01734 CALL hmpdef(13,-prange,prange, & ! book 01735 'Pulls of single (global) measurement') 01736 CALL hmpdef(14,-prange,prange, & ! book 01737 'Pulls of single (local) measurement') 01738 CALL hmpdef(4,0.0,0.0,'Chi^2/Ndf after local fit') 01739 CALL gmpdef(4,5,'location, dispersion (res.) vs record nr') 01740 CALL gmpdef(5,5,'location, dispersion (pull) vs record nr') 01741 01742 ! WRITE(*,*) 'LOOPN ', NLOOPN, ' executing ICALCM=', ICALCM 01743 01744 ! reset 01745 01746 globalVector=0.0D0 ! reset rhs vector IGVEC 01747 IF(icalcm == 1) THEN 01748 globalMatD=0.0D0 01749 globalMatF=0. 01750 IF (metsol >= 3) matPreCond=0.0D0 01751 END IF 01752 01753 IF(nloopn == 2) CALL hmpdef(6,0.0,0.0,'Down-weight fraction') 01754 01755 newite=.FALSE. 01756 IF(iterat /= lastit) THEN ! new iteration 01757 newite=.TRUE. 01758 funref=fvalue 01759 IF(nloopn > 1) THEN 01760 nrej=nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) 01761 ! CALL MEND 01762 IF(iterat == 1) THEN 01763 chicut=chirem 01764 ELSE IF(iterat >= 1) THEN 01765 chicut=SQRT(chicut) 01766 IF(chicut /= 0.0.AND.chicut < 1.5) chicut=1.0 01767 IF(chicut /= 0.0.AND.nrej == 0) chicut=1.0 01768 END IF 01769 END IF 01770 ! WRITE(*,111) ! header line 01771 END IF 01772 01773 DO i=0,3 01774 nrejec(i)=0 ! reset reject counter 01775 END DO 01776 DO k=3,6 01777 writeBufferHeader(k)=0 ! cache usage 01778 writeBufferHeader(-k)=0 01779 END DO 01780 ! statistics per binary file 01781 DO i=1,nfilb 01782 jfd(i)=0 01783 cfd(i)=0.0 01784 dfd(i)=0 01785 END DO 01786 ! ----- read next data ---------------------------------------------- 01787 DO 01788 CALL peread(nr) ! read records 01789 CALL peprep(1) ! prepare records 01790 ndfs =0 01791 sndf =0.0D0 01792 dchi2s=0.0D0 01793 CALL loopbf(nrejec,ndfs,sndf,dchi2s,nfiles,jfd,cfd,dfd) 01794 ndfsum=ndfsum+ndfs 01795 sumndf=sumndf+sndf 01796 CALL addsum(dchi2s) 01797 IF (nr <= 0) EXIT ! next block of events ? 01798 END DO 01799 ! sum up RHS (over threads) once (reduction in LOOPBF: summation for each block) 01800 ioffb=0 01801 DO ipr=2,mthrd 01802 ioffb=ioffb+lenGlobalVec 01803 DO k=1,lenGlobalVec 01804 globalVector(k)=globalVector(k)+globalVector(ioffb+k) 01805 END DO 01806 END DO 01807 01808 IF (icalcm == 1) THEN 01809 ! PRINT *, ' cache/w ',(writeBufferHeader(-K),K=3,6),(writeBufferHeader(K),K=3,6) 01810 nparl=writeBufferHeader(3) 01811 ncrit=writeBufferHeader(4) 01812 used=FLOAT(writeBufferHeader(-5))/FLOAT(writeBufferHeader(-3))*0.1 01813 usei=FLOAT(writeBufferHeader(5))/FLOAT(writeBufferHeader(3))*0.1 01814 peakd=FLOAT(writeBufferHeader(-6))*0.1 01815 peaki=FLOAT(writeBufferHeader(6))*0.1 01816 WRITE(*,111) nparl,ncrit,usei,used,peaki,peakd 01817 111 FORMAT(' Write cache usage (#flush,#overrun,<levels>,', & 01818 'peak(levels))'/2I7,',',4(f6.1,'%')) 01819 END IF 01820 01821 ! ----- after end-of-data add contributions from pre-sigma --------- 01822 01823 IF(nloopn == 1) THEN 01824 ! plot diagonal elements 01825 elmt=0.0 01826 DO i=1,nvgb ! diagonal elements 01827 ii=0 01828 IF(matsto == 1) THEN 01829 ii=i 01830 ii=(ii*ii+ii)/2 01831 END IF 01832 IF(matsto == 2) ii=i 01833 IF(matsto == 3) ii=i 01834 IF(ii /= 0) THEN 01835 elmt=SNGL(globalMatD(ii)) 01836 IF(elmt > 0.0) CALL hmpent(23,1.0/SQRT(elmt)) 01837 END IF 01838 END DO 01839 END IF 01840 01841 01842 01843 ! add pre-sigma contributions to matrix diagonal 01844 01845 ! WRITE(*,*) 'Adding to diagonal ICALCM IND6',ICALCM,IND6 01846 01847 IF(icalcm == 1) THEN 01848 DO ivgb=1,nvgb ! add evtl. pre-sigma 01849 ! WRITE(*,*) 'Index ',IVGB,IVGB,QM(IND6+IVGB) 01850 IF(globalParPreWeight(ivgb) /= 0.0) THEN 01851 IF(ivgb > 0) CALL mupdat(ivgb,ivgb,DBLE(globalParPreWeight(ivgb))) 01852 END IF 01853 END DO 01854 END IF 01855 01856 CALL hmpwrt(23) 01857 CALL hmpwrt(24) 01858 CALL hmpwrt(25) 01859 CALL hmpwrt(26) 01860 01861 01862 ! add regularization term to F and to rhs -------------------------- 01863 01864 ! WRITE(*,*) 'NREGUL ',NREGUL,NLOOPN 01865 01866 IF(nregul /= 0) THEN ! add regularization term to F and to rhs 01867 DO ivgb=1,nvgb 01868 itgbi=globalParVarToTotal(ivgb) ! global parameter index 01869 globalVector(ivgb)=globalVector(ivgb) -globalParameter(itgbi)*globalParPreWeight(ivgb) 01870 adder=globalParPreWeight(ivgb)*globalParameter(itgbi)**2 01871 CALL addsum(adder) 01872 END DO 01873 END IF 01874 01875 01876 ! ----- add contributions from "measurement" ----------------------- 01877 01878 01879 i=1 01880 DO WHILE (i <= lenMeasurements) 01881 rhs=listMeasurements(i )%value ! right hand side 01882 sgm=listMeasurements(i+1)%value ! sigma parameter 01883 i=i+2 01884 weight=0.0 01885 IF(sgm > 0.0) weight=1.0/sgm**2 01886 01887 dsum=-rhs 01888 01889 ! loop over label/factor pairs 01890 ia=i 01891 DO 01892 i=i+1 01893 IF(i > lenMeasurements) EXIT 01894 IF(listMeasurements(i)%label == 0) EXIT 01895 END DO 01896 ib=i-1 01897 01898 DO j=ia,ib 01899 factrj=listMeasurements(j)%value 01900 itgbij=inone(listMeasurements(j)%label) ! total parameter index 01901 IF(itgbij /= 0) THEN 01902 dsum=dsum+factrj*SNGL(globalParameter(itgbij)) ! residuum 01903 END IF 01904 ! add to vector 01905 ivgbij=0 01906 IF(itgbij /= 0) ivgbij=globalParLabelIndex(2,itgbij) ! variable-parameter index 01907 IF(ivgbij > 0) THEN 01908 globalVector(ivgbij)=globalVector(ivgbij) -weight*dsum*factrj ! vector 01909 END IF 01910 01911 IF(icalcm == 1.AND.ivgbij > 0) THEN 01912 DO k=ia,j 01913 factrk=listMeasurements(k)%value 01914 itgbik=inone(listMeasurements(k)%label) ! total parameter index 01915 ! add to matrix 01916 ivgbik=0 01917 IF(itgbik /= 0) ivgbik=globalParLabelIndex(2,itgbik) ! variable-parameter index 01918 IF(ivgbij > 0.AND.ivgbik > 0) THEN ! 01919 CALL mupdat(ivgbij,ivgbik,DBLE(weight*factrj*factrk)) 01920 END IF 01921 END DO 01922 END IF 01923 END DO 01924 01925 adder=weight*dsum**2 01926 CALL addsum(adder) ! accumulate chi-square 01927 01928 END DO 01929 01930 ! ----- printout --------------------------------------------------- 01931 01932 01933 CALL getsum(fvalue) ! get accurate sum (Chi^2) 01934 flines=0.5D0*fvalue ! Likelihood function value 01935 rloop=iterat+0.01*nloopn 01936 actfun=SNGL(funref-fvalue) 01937 IF(nloopn == 1) actfun=0.0 01938 ngras=nint(angras) 01939 ratae=0.0 !!! 01940 IF(delfun /= 0.0) THEN 01941 ratae=MIN(99.9,actfun/delfun) !!! 01942 ratae=MAX(-99.9,ratae) 01943 END IF 01944 01945 ! rejects ... 01946 01947 nrej =nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) 01948 IF(nloopn == 1) THEN 01949 IF(nrej /= 0) THEN 01950 WRITE(*,*) ' ' 01951 WRITE(*,*) 'Data rejected in initial loop:' 01952 WRITE(*,*) ' ', & 01953 nrejec(0), ' (rank deficit/NaN) ',nrejec(1),' (Ndf=0) ', & 01954 nrejec(2), ' (huge) ',nrejec(3),' (large)' 01955 END IF 01956 END IF 01957 ! IF(NREJEC(1)+NREJEC(2)+NREJEC(3).NE.0) THEN 01958 ! WRITE(LUNLOG,*) 'Data rejected in initial loop:',NREJEC(1), 01959 ! + ' (Ndf=0) ',NREJEC(2),' (huge) ',NREJEC(3),' (large)' 01960 ! END IF 01961 01962 01963 IF(newite.AND.iterat == 2) THEN 01964 IF(nrecpr /= 0.OR.nrecp2 /= 0) nrecer=nrec3 01965 IF(nrecpr < 0) THEN 01966 nrecpr=nrec1 01967 END IF 01968 IF(nrecp2 < 0) THEN 01969 nrecp2=nrec2 01970 END IF 01971 END IF 01972 01973 IF(nloopn <= 2) THEN 01974 IF(nhistp /= 0) THEN 01975 ! CALL HMPRNT(3) ! scaled residual of single measurement 01976 ! CALL HMPRNT(12) ! scaled residual of single measurement 01977 ! CALL HMPRNT(4) ! chi^2/Ndf 01978 END IF 01979 CALL hmpwrt(3) 01980 CALL hmpwrt(12) 01981 CALL hmpwrt(4) 01982 CALL gmpwrt(4) ! location, dispersion (res.) as a function of record nr 01983 IF (nloopn <= lfitnp) THEN 01984 CALL hmpwrt(13) 01985 CALL hmpwrt(14) 01986 CALL gmpwrt(5) ! location, dispersion (pull) as a function of record nr 01987 END IF 01988 END IF 01989 ! IF(NLOOPN.EQ.2.AND.NHISTP.NE.0) CALL HMPRNT(6) 01990 IF(nloopn == 2) CALL hmpwrt(6) 01991 IF(nloopn <= 1) THEN 01992 ! IF(NHISTP.NE.0) CALL HMPRNT(5) ! number of degrees of freedom 01993 ! IF(NHISTP.NE.0) CALL HMPRNT(11) ! Nlocal 01994 CALL hmpwrt(5) 01995 CALL hmpwrt(11) 01996 END IF 01997 01998 ! local fit: band matrix structure !? 01999 IF (nloopn == 1.AND.nbndr > 0) THEN 02000 DO lun=6,8,2 02001 WRITE(lun,*) ' ' 02002 WRITE(lun,*) ' === local fits have bordered band matrix structure ===' 02003 WRITE(lun,101) ' NBNDR',nbndr,'number of records' 02004 WRITE(lun,101) ' NBDRX',nbdrx,'max border size' 02005 WRITE(lun,101) ' NBNDX',nbndx,'max band width' 02006 END DO 02007 END IF 02008 02009 lastit=iterat 02010 02011 101 FORMAT(1X,a6,' =',i10,' = ',a) 02012 ! 101 FORMAT(' LOOPN',I6,' Function value',F22.8,10X,I6,' records') 02013 ! 102 FORMAT(' incl. constraint penalty',F22.8) 02014 ! 103 FORMAT(I13,3X,A,G12.4) 02015 END SUBROUTINE loopn ! loop with fits 02016 02020 02021 SUBROUTINE ploopa(lunp) 02022 USE mpmod 02023 02024 IMPLICIT NONE 02025 02026 INTEGER, INTENT(IN) :: lunp 02027 ! .. 02028 WRITE(lunp,*) ' ' 02029 WRITE(lunp,101) ! header line 02030 WRITE(lunp,102) ! header line 02031 101 FORMAT(' it fc',' fcn_value dfcn_exp slpr costh iit st', & 02032 ' ls step cutf',1X,'rejects hmmsec FMS') 02033 102 FORMAT(' -- --',' ----------- -------- ---- ----- --- --', & 02034 ' -- ---- ----',1X,'------- ------ ---') 02035 RETURN 02036 END SUBROUTINE ploopa ! title for iteration 02037 02041 02042 SUBROUTINE ploopb(lunp) 02043 USE mpmod 02044 02045 IMPLICIT NONE 02046 INTEGER :: ma 02047 INTEGER :: minut 02048 INTEGER :: nfa 02049 INTEGER :: nhour 02050 INTEGER :: nrej 02051 INTEGER :: nsecnd 02052 REAL :: ratae 02053 REAL :: rstb 02054 REAL :: secnd 02055 REAL :: slopes(3) 02056 REAL :: steps(3) 02057 REAL, DIMENSION(2) :: ta 02058 02059 INTEGER, INTENT(IN) :: lunp 02060 02061 CHARACTER (LEN=4):: ccalcm(4) 02062 DATA ccalcm / ' end',' S', ' F ',' FMS' / 02063 SAVE 02064 02065 nrej=nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) ! rejects 02066 IF(nrej > 9999999) nrej=9999999 02067 CALL etime(ta,rstb) 02068 deltim=rstb-rstart 02069 CALL sechms(deltim,nhour,minut,secnd) ! time 02070 nsecnd=nint(secnd) 02071 IF(iterat == 0) THEN 02072 WRITE(lunp,103) iterat,nloopn,fvalue, & 02073 chicut,nrej,nhour,minut,nsecnd,ccalcm(lcalcm) 02074 ELSE 02075 CALL ptlopt(nfa,ma,slopes,steps) ! slopes steps 02076 ratae=ABS(slopes(2)/slopes(1)) 02077 stepl=steps(2) 02078 WRITE(lunp,104) iterat,nloopn,fvalue,delfun,ratae,angras, & 02079 iitera,istopa,lsinfo,stepl, chicut,nrej,nhour,minut,nsecnd,ccalcm(lcalcm) 02080 END IF 02081 103 FORMAT(i3,i3,e12.5,38X,f5.1, 1X,i7, i2,i2,i3,a4) 02082 104 FORMAT(i3,i3,e12.5,1X,e8.2,f6.3,f6.3,i5,2I3,f6.2,f5.1, & 02083 1X,i7, i2,i2,i3,a4) 02084 RETURN 02085 END SUBROUTINE ploopb ! iteration line 02086 02090 02091 SUBROUTINE ploopc(lunp) 02092 USE mpmod 02093 02094 IMPLICIT NONE 02095 INTEGER :: ma 02096 INTEGER :: minut 02097 INTEGER :: nfa 02098 INTEGER :: nhour 02099 INTEGER :: nrej 02100 INTEGER :: nsecnd 02101 REAL :: ratae 02102 REAL :: rstb 02103 REAL :: secnd 02104 REAL :: slopes(3) 02105 REAL :: steps(3) 02106 REAL, DIMENSION(2) :: ta 02107 02108 INTEGER, INTENT(IN) :: lunp 02109 CHARACTER (LEN=4):: ccalcm(4) 02110 DATA ccalcm / ' end',' S', ' F ',' FMS' / 02111 SAVE 02112 02113 nrej=nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) ! rejects 02114 IF(nrej > 9999999) nrej=9999999 02115 CALL etime(ta,rstb) 02116 deltim=rstb-rstart 02117 CALL sechms(deltim,nhour,minut,secnd) ! time 02118 nsecnd=nint(secnd) 02119 CALL ptlopt(nfa,ma,slopes,steps) ! slopes steps 02120 ratae=ABS(slopes(2)/slopes(1)) 02121 stepl=steps(2) 02122 WRITE(lunp,105) nloopn,fvalue, ratae,lsinfo, & 02123 stepl,nrej,nhour,minut,nsecnd,ccalcm(lcalcm) 02124 105 FORMAT(3X,i3,e12.5,9X, f6.3,14X,i3,f6.2,6X, i7, i2,i2,i3,a4) 02125 RETURN 02126 02127 END SUBROUTINE ploopc ! sub-iteration line 02128 02132 02133 SUBROUTINE ploopd(lunp) 02134 USE mpmod 02135 IMPLICIT NONE 02136 INTEGER :: minut 02137 INTEGER :: nhour 02138 INTEGER :: nsecnd 02139 REAL :: rstb 02140 REAL :: secnd 02141 REAL, DIMENSION(2) :: ta 02142 02143 INTEGER, INTENT(IN) :: lunp 02144 CHARACTER (LEN=4):: ccalcm(4) 02145 DATA ccalcm / ' end',' S', ' F ',' FMS' / 02146 SAVE 02147 CALL etime(ta,rstb) 02148 deltim=rstb-rstart 02149 CALL sechms(deltim,nhour,minut,secnd) ! time 02150 nsecnd=IFIX(secnd+0.5) 02151 02152 WRITE(lunp,106) nhour,minut,nsecnd,ccalcm(lcalcm) 02153 106 FORMAT(69X,i2,i2,i3,a4) 02154 RETURN 02155 END SUBROUTINE ploopd 02156 02158 SUBROUTINE explfc(lunit) 02159 IMPLICIT NONE 02160 INTEGER:: lunit 02161 WRITE(lunit,*) ' ' 02162 WRITE(lunit,102) 'Explanation of iteration table' 02163 WRITE(lunit,102) '==============================' 02164 WRITE(lunit,101) 'it', & 02165 'iteration number. Global parameters are improved for it > 0.' 02166 WRITE(lunit,102) 'First function evaluation is called iteraton 0.' 02167 WRITE(lunit,101) 'fc', 'number of function evaluations.' 02168 WRITE(lunit,101) 'fcn_value', 'value of 2 x Likelihood function (LF).' 02169 WRITE(lunit,102) 'The final value is the chi^2 value of the fit and should' 02170 WRITE(lunit,102) 'be about equal to the NDF (see below).' 02171 WRITE(lunit,101) 'dfcn_exp', & 02172 'expected reduction of the value of the Likelihood function (LF)' 02173 WRITE(lunit,101) 'slpr', 'ratio of the actual slope to inital slope.' 02174 WRITE(lunit,101) 'costh', & 02175 'cosine of angle between search direction and -gradient' 02176 WRITE(lunit,101) 'iit', & 02177 'number of internal iterations in GMRES/MINRES algorithmus' 02178 WRITE(lunit,101) 'st', 'stop code of GMRES/MINRES algorithmus' 02179 WRITE(lunit,102) '< 0: rhs is very special, with beta2 = 0' 02180 WRITE(lunit,102) '= 0: rhs b = 0, i.e. the exact solution is x = 0' 02181 WRITE(lunit,102) '= 1 requested accuracy achieved, as determined by rtol' 02182 WRITE(lunit,102) '= 2 reasonable accuracy achieved, given eps' 02183 WRITE(lunit,102) '= 3 x has converged to an eigenvector' 02184 WRITE(lunit,102) '= 4 matrix ill-conditioned (Acond has exceeded 0.1/eps)' 02185 WRITE(lunit,102) '= 5 the iteration limit was reached' 02186 WRITE(lunit,102) '= 6 Matrix x vector does not define a symmetric matrix' 02187 WRITE(lunit,102) '= 7 Preconditioner does not define a symmetric matrix' 02188 WRITE(lunit,101) 'ls', 'line search info' 02189 WRITE(lunit,102) '< 0 recalculate function' 02190 WRITE(lunit,102) '= 0: N or STP lt 0 or step not descending' 02191 WRITE(lunit,102) '= 1: Linesearch convergence conditions reached' 02192 WRITE(lunit,102) '= 2: interval of uncertainty at lower limit' 02193 WRITE(lunit,102) '= 3: max nr of line search calls reached' 02194 WRITE(lunit,102) '= 4: step at the lower bound' 02195 WRITE(lunit,102) '= 5: step at the upper bound' 02196 WRITE(lunit,102) '= 6: rounding error limitation' 02197 WRITE(lunit,101) 'step', & 02198 'the factor for the Newton step during the line search. Usually' 02199 WRITE(lunit,102) & 02200 'a value of 1 gives a sufficient reduction of the LF. Oherwise' 02201 WRITE(lunit,102) 'other step values are tried.' 02202 WRITE(lunit,101) 'cutf', & 02203 'cut factor. Local fits are rejected, if their chi^2 value' 02204 WRITE(lunit,102) & 02205 'is larger than the 3-sigma chi^2 value times the cut factor.' 02206 WRITE(lunit,102) 'A cut factor of 1 is used finally, but initially a larger' 02207 WRITE(lunit,102) 'factor may be used. A value of 0.0 means no cut.' 02208 WRITE(lunit,101) 'rejects', 'total number of rejected local fits.' 02209 WRITE(lunit,101) 'hmmsec', 'the time in hours (h), minutes (mm) and seconds.' 02210 WRITE(lunit,101) 'FMS', 'calculation of Function value, Matrix, Solution.' 02211 WRITE(lunit,*) ' ' 02212 02213 101 FORMAT(a9,' = ',a) 02214 102 FORMAT(13X,a) 02215 END SUBROUTINE explfc 02216 02224 02225 SUBROUTINE mupdat(i,j,add) ! 02226 USE mpmod 02227 02228 IMPLICIT NONE 02229 02230 INTEGER, INTENT(IN) :: i 02231 INTEGER, INTENT(IN) :: j 02232 DOUBLE PRECISION, INTENT(IN) :: add 02233 02234 INTEGER(KIND=large):: ijadd 02235 INTEGER(KIND=large):: ia 02236 INTEGER(KIND=large):: ja 02237 INTEGER(KIND=large):: ij 02238 ! ... 02239 IF(i <= 0.OR.j <= 0) RETURN 02240 ia=MAX(i,j) ! larger 02241 ja=MIN(i,j) ! smaller 02242 ij=0 02243 IF(matsto == 1) THEN ! full symmetric matrix 02244 ij=ja+(ia*ia-ia)/2 ! ISYM index 02245 globalMatD(ij)=globalMatD(ij)+add 02246 ELSE IF(matsto == 2) THEN ! sparse symmetric matrix 02247 ij=ijadd(i,j) ! inline code requires same time 02248 IF (ij == 0) RETURN ! pair is suppressed 02249 IF (ij > 0) THEN 02250 globalMatD(ij)=globalMatD(ij)+add 02251 ELSE 02252 globalMatF(-ij)=globalMatF(-ij)+SNGL(add) 02253 END IF 02254 END IF 02255 IF(metsol >= 3) THEN 02256 IF(mbandw > 0) THEN ! for Cholesky decomposition 02257 IF(ia <= nvgb) THEN ! variable global parameter 02258 ij=indPreCond(ia)-ia+ja 02259 IF(ia > 1.AND.ij <= indPreCond(ia-1)) ij=0 02260 IF(ij /= 0) matPreCond(ij)=matPreCond(ij)+add 02261 IF(ij < 0.OR.ij > size(matPreCond)) STOP 'mupdat: bad index' 02262 ELSE ! Lagrange multiplier 02263 ij=indPreCond(nvgb)+(ia-nvgb-1)*nvgb+ja 02264 IF(ij /= 0) matPreCond(ij)=matPreCond(ij)+add 02265 END IF 02266 ELSE IF(mbandw == 0) THEN ! default preconditioner 02267 IF(ia <= nvgb) THEN ! variable global parameter 02268 IF(ja == ia) matPreCond(ia)=matPreCond(ia)+add ! diag 02269 ELSE ! Lagrange multiplier 02270 ij=nvgb+(ia-nvgb-1)*nvgb+ja 02271 IF(ij /= 0) matPreCond(ij)=matPreCond(ij)+add 02272 END IF 02273 END IF 02274 END IF 02275 END SUBROUTINE mupdat 02276 02277 02307 02308 SUBROUTINE loopbf(nrej,ndfs,sndf,dchi2s, numfil,naccf,chi2f,ndff) 02309 USE mpmod 02310 02311 IMPLICIT NONE 02312 REAL :: cauchy 02313 REAL :: chichi 02314 REAL :: chlimt 02315 REAL :: chndf 02316 REAL :: chuber 02317 REAL :: down 02318 REAL :: glder 02319 REAL :: pull 02320 REAL :: r1 02321 REAL :: r2 02322 REAL :: rec 02323 REAL :: rerr 02324 REAL :: resid 02325 REAL :: resing 02326 REAL :: resmax 02327 REAL :: rmeas 02328 REAL :: rmloc 02329 REAL :: suwt 02330 REAL :: used 02331 REAL :: wght 02332 REAL :: wrc 02333 REAL :: chindl 02334 INTEGER :: i 02335 INTEGER :: ia 02336 INTEGER :: ib 02337 INTEGER :: ibuf 02338 INTEGER :: ichunk 02339 INTEGER :: icmn 02340 INTEGER :: icost 02341 INTEGER :: id 02342 INTEGER :: idiag 02343 INTEGER :: iext 02344 INTEGER :: ij 02345 INTEGER :: ije 02346 INTEGER :: ijn 02347 INTEGER :: ijsym 02348 INTEGER :: ik 02349 INTEGER :: ike 02350 INTEGER :: im 02351 INTEGER :: imeas 02352 INTEGER :: in 02353 INTEGER :: inder 02354 INTEGER :: inv 02355 INTEGER :: ioffb 02356 INTEGER :: ioffc 02357 INTEGER :: ioffd 02358 INTEGER :: ioffe 02359 INTEGER :: ioffi 02360 INTEGER :: iprdbg 02361 INTEGER :: iproc 02362 INTEGER :: irow 02363 INTEGER :: isfrst 02364 INTEGER :: islast 02365 INTEGER :: ist 02366 INTEGER :: iter 02367 INTEGER :: itgbi 02368 INTEGER :: ivgbj 02369 INTEGER :: ivgbk 02370 INTEGER :: j 02371 INTEGER :: ja 02372 INTEGER :: jb 02373 INTEGER :: jk 02374 INTEGER :: jn 02375 INTEGER :: joffd 02376 INTEGER :: joffi 02377 INTEGER :: jproc 02378 INTEGER :: jsp 02379 INTEGER :: k 02380 INTEGER :: kbdr 02381 INTEGER :: kbdrx 02382 INTEGER :: kbnd 02383 INTEGER :: kfl 02384 INTEGER :: kx 02385 INTEGER :: mbdr 02386 INTEGER :: mbnd 02387 INTEGER :: nalc 02388 INTEGER :: nalg 02389 INTEGER :: nan 02390 INTEGER :: nb 02391 INTEGER :: ndf 02392 INTEGER :: ndfsd 02393 INTEGER :: ndown 02394 INTEGER :: neq 02395 INTEGER :: nfred 02396 INTEGER :: nfrei 02397 INTEGER :: ngg 02398 INTEGER :: nprdbg 02399 INTEGER :: nrank 02400 INTEGER :: nrc 02401 INTEGER :: nst 02402 INTEGER :: nter 02403 INTEGER :: nweig 02404 02405 INTEGER, INTENT(IN OUT) :: nrej(0:3) 02406 INTEGER, INTENT(IN OUT) :: ndfs 02407 DOUBLE PRECISION, INTENT(IN OUT) :: sndf 02408 DOUBLE PRECISION, INTENT(IN OUT) :: dchi2s 02409 INTEGER, INTENT(IN) :: numfil 02410 INTEGER, INTENT(IN OUT) :: naccf(numfil) 02411 REAL, INTENT(IN OUT) :: chi2f(numfil) 02412 INTEGER, INTENT(IN OUT) :: ndff(numfil) 02413 02414 DOUBLE PRECISION::dch2sd 02415 DOUBLE PRECISION:: dchi2 02416 DOUBLE PRECISION::dvar 02417 DOUBLE PRECISION:: dw1 02418 DOUBLE PRECISION::dw2 02419 DOUBLE PRECISION::summ 02420 02421 !$ INTEGER OMP_GET_THREAD_NUM 02422 02423 LOGICAL:: lprnt 02424 LOGICAL::lhist 02425 CHARACTER (LEN=3):: chast 02426 DATA chuber/1.345/ ! constant for Huber down-weighting 02427 DATA cauchy/2.3849/ ! constant for Cauchy down-weighting 02428 SAVE chuber,cauchy 02429 ! ... 02430 ijsym(i,j)=MIN(i,j)+(MAX(i,j)*MAX(i,j)-MAX(i,j))/2 02431 isfrst(ibuf)=readBufferPointer(ibuf)+1 02432 islast(ibuf)=readBufferDataI(readBufferPointer(ibuf)) 02433 inder(i)=readBufferDataI(i) 02434 glder(i)=readBufferDataF(i) 02435 02436 ichunk=MIN((numReadBuffer+mthrd-1)/mthrd/32+1,256) 02437 ! reset header, 3 words per thread: 02438 ! number of entries, offset to data, indices 02439 writeBufferInfo=0 02440 writeBufferData=0. 02441 ! IF (ICALCM.EQ.1) print *, ' ICHUNK ', MQ(IDOT3),MTHRD,ICHUNK 02442 nprdbg=0 02443 iprdbg=-1 02444 02445 ! parallelize record loop 02446 ! private copy of NDFS,.. for each thread, combined at end, init with 0. 02447 !$OMP PARALLEL DO & 02448 !$OMP DEFAULT(PRIVATE) & 02449 !$OMP SHARED(numReadBuffer,readBufferPointer,readBufferDataI, & 02450 !$OMP readBufferDataF,writeBufferHeader,writeBufferInfo, & 02451 !$OMP writeBufferData,writeBufferIndices,writeBufferUpdates,globalVector, & 02452 !$OMP globalParameter,globalParLabelIndex,globalIndexUsage,backIndexUsage, & 02453 !$OMP NAGB,NVGB,NAGBN,ICALCM,ICHUNK,NLOOPN,NRECER,NPRDBG,IPRDBG, & 02454 !$OMP NEWITE,CHICUT,LHUBER,CHUBER,ITERAT,NRECPR,MTHRD, & 02455 !$OMP DWCUT,CHHUGE,NRECP2,CAUCHY,LFITNP,LFITBB) & 02456 !$OMP REDUCTION(+:NDFS,SNDF,DCHI2S,NREJ,NBNDR,NACCF,CHI2F,NDFF) & 02457 !$OMP REDUCTION(MAX:NBNDX,NBDRX) & 02458 !$OMP REDUCTION(MIN:NREC3) & 02459 !$OMP SCHEDULE(DYNAMIC,ICHUNK) 02460 DO ibuf=1,numReadBuffer ! buffer for current record 02461 nrc=readBufferDataI(isfrst(ibuf)-2) ! record 02462 kfl=NINT(readBufferDataF(isfrst(ibuf)-1)) ! file 02463 wrc=readBufferDataF(isfrst(ibuf)-2) ! weight 02464 dw1=DBLE(wrc) 02465 dw2=DSQRT(dw1) 02466 02467 iproc=0 02468 !$ IPROC=OMP_GET_THREAD_NUM() ! thread number 02469 ioffb=nagb*iproc ! offset 'f'. 02470 ioffc=nagbn*iproc ! offset 'c'. 02471 ioffe=nvgb*iproc ! offset 'e' 02472 ioffd=writeBufferHeader(-1)*iproc+writeBufferInfo(2,iproc+1) ! offset data 02473 ioffi=writeBufferHeader(1)*iproc+writeBufferInfo(3,iproc+1)+2 ! offset indices 02474 ! ----- reset ------------------------------------------------------ 02475 lprnt=.FALSE. 02476 lhist=(iproc == 0) 02477 REC=nrc ! floating point value 02478 IF(nloopn == 1.AND.MOD(nrc,100000) == 0) THEN 02479 WRITE(*,*) 'Record',nrc,' ... still reading' 02480 END IF 02481 02482 ! printout/debug only for one thread at a time 02483 02484 02485 ! flag for record printout ----------------------------------------- 02486 02487 lprnt=.FALSE. 02488 IF(newite.AND.(iterat == 1.OR.iterat == 3)) THEN 02489 IF(nrc == nrecpr) lprnt=.TRUE. 02490 IF(nrc == nrecp2) lprnt=.TRUE. 02491 IF(nrc == nrecer) lprnt=.TRUE. 02492 END IF 02493 IF (lprnt)THEN 02494 !$OMP ATOMIC 02495 nprdbg=nprdbg+1 ! number of threads with debug 02496 IF (nprdbg == 1) iprdbg=iproc ! first thread with debug 02497 IF (iproc /= iprdbg) lprnt=.FALSE. 02498 ! print *, ' LPRNT ', NRC, NPRDBG, IPRDBG, IPROC, LPRNT 02499 END IF 02500 IF(lprnt) THEN 02501 WRITE(1,*) ' ' 02502 WRITE(1,*) '------------------ Loop',nloopn, & 02503 ': Printout for record',nrc,iproc 02504 WRITE(1,*) ' ' 02505 END IF 02506 02507 ! ----- print data ------------------------------------------------- 02508 02509 IF(lprnt) THEN 02510 imeas=0 ! local derivatives 02511 ist=isfrst(ibuf) 02512 nst=islast(ibuf) 02513 DO ! loop over measurements 02514 CALL isjajb(nst,ist,ja,jb,jsp) 02515 IF(ja == 0) EXIT 02516 IF(imeas == 0) WRITE(1,1121) 02517 imeas=imeas+1 02518 WRITE(1,1122) imeas,glder(ja),glder(jb), & 02519 (inder(ja+j),glder(ja+j),j=1,jb-ja-1) 02520 END DO 02521 1121 FORMAT(/'Measured value and local derivatives'/ & 02522 ' i measured std_dev index...derivative ...') 02523 1122 FORMAT(i3,2G12.4,3(i3,g12.4)/(27X,3(i3,g12.4))) 02524 02525 imeas=0 ! global derivatives 02526 ist=isfrst(ibuf) 02527 nst=islast(ibuf) 02528 DO ! loop over measurements 02529 CALL isjajb(nst,ist,ja,jb,jsp) 02530 IF(ja == 0) EXIT 02531 IF(imeas == 0) WRITE(1,1123) 02532 imeas=imeas+1 02533 IF (jb < ist) THEN 02534 IF(ist-jb > 2) THEN 02535 WRITE(1,1124) imeas,(globalParLabelIndex(1,inder(jb+j)),inder(jb+j), & 02536 globalParLabelIndex(2,inder(jb+j)),glder(jb+j),j=1,ist-jb) 02537 ELSE 02538 WRITE(1,1125) imeas,(globalParLabelIndex(1,inder(jb+j)),inder(jb+j), & 02539 globalParLabelIndex(2,inder(jb+j)),glder(jb+j),j=1,ist-jb) 02540 END IF 02541 END IF 02542 END DO 02543 1123 FORMAT(/'Global derivatives'/ & 02544 ' i label gindex vindex derivative ...') 02545 1124 FORMAT(i3,2(i9,i7,i7,g12.4)/(3X,2(i9,i7,i7,g12.4))) 02546 1125 FORMAT(i3,2(i9,i7,i7,g12.4)) 02547 END IF 02548 02549 ! ----- first loop ------------------------------------------------- 02550 ! ------ prepare local fit ------ 02551 ! count local and global derivates 02552 ! subtract actual alignment parameters from the measured data 02553 02554 IF(lprnt) THEN 02555 WRITE(1,*) ' ' 02556 WRITE(1,*) 'Data corrections using values of global parameters' 02557 WRITE(1,*) '==================================================' 02558 WRITE(1,101) 02559 END IF 02560 nalg=0 ! count number of global derivatives 02561 nalc=0 ! count number of local derivatives 02562 imeas=0 02563 ist=isfrst(ibuf) 02564 nst=islast(ibuf) 02565 DO ! loop over measurements 02566 CALL isjajb(nst,ist,ja,jb,jsp) 02567 IF(ja == 0) EXIT 02568 rmeas=glder(ja) ! data 02569 ! subtract global ... from measured value 02570 DO j=1,ist-jb ! global parameter loop 02571 itgbi=inder(jb+j) ! global parameter label 02572 rmeas=rmeas-glder(jb+j)*SNGL(globalParameter(itgbi)) ! subtract !!! reversed 02573 IF (icalcm == 1) THEN 02574 ij=globalParLabelIndex(2,itgbi) ! index of variable global parameter 02575 IF(ij > 0) THEN 02576 ijn=backIndexUsage(ioffe+ij) ! get index of index 02577 IF(ijn == 0) THEN ! not yet included 02578 nalg=nalg+1 ! count 02579 globalIndexUsage(ioffc+nalg)=ij ! store global index 02580 backIndexUsage(ioffe+ij)=nalg ! store back index 02581 END IF 02582 END IF 02583 END IF 02584 END DO 02585 IF(lprnt) THEN 02586 imeas=imeas+1 02587 IF (jb < ist) WRITE(1,102) imeas,glder(ja),rmeas,glder(jb) 02588 END IF 02589 readBufferDataF(ja)=rmeas ! global contribution subtracted 02590 DO j=1,jb-ja-1 ! local parameter loop 02591 ij=inder(ja+j) 02592 nalc=MAX(nalc,ij) ! number of local parameters 02593 END DO 02594 END DO 02595 101 FORMAT(' index measvalue corrvalue sigma') 02596 102 FORMAT(i6,2X,2G12.4,' +-',g12.4) 02597 02598 IF(nalc <= 0) GO TO 90 02599 02600 ngg=(nalg*nalg+nalg)/2 02601 IF (icalcm == 1) THEN 02602 DO k=1,nalg*nalc 02603 localGlobalMatrix(k)=0.0D0 ! reset global-local matrix 02604 END DO 02605 writeBufferIndices(ioffi-1)=nrc ! index header: 02606 writeBufferIndices(ioffi )=nalg ! event number, number of global par 02607 CALL sort1k(globalIndexUsage(ioffc+1),nalg) ! sort global par. 02608 DO k=1,nalg 02609 iext=globalIndexUsage(ioffc+k) 02610 writeBufferIndices(ioffi+k)=iext ! global par indices 02611 backIndexUsage(ioffe+iext)=k ! update back index 02612 END DO 02613 DO k=1,ngg 02614 writeBufferUpdates(ioffd+k)=0.0D0 ! reset global-global matrix 02615 END DO 02616 END IF 02617 ! ----- iteration start and check --------------------------------- 02618 02619 nter=1 ! first loop without down-weighting 02620 IF(nloopn /= 1.AND.lhuber /= 0) nter=lhuber 02621 02622 ! check matrix for bordered band structure (MBDR+MBND+1 <= NALC) 02623 mbnd=-1 02624 mbdr=nalc 02625 DO i=1, nalc 02626 ibandh(i)=0 02627 END DO 02628 irow=1 02629 idiag=1 02630 ndfsd=0 02631 02632 iter=0 02633 resmax=0.0 02634 DO WHILE(iter < nter) ! outlier suppresssion iteration loop 02635 iter=iter+1 02636 resmax=0.0 02637 IF(lprnt) THEN 02638 WRITE(1,*) ' ' 02639 WRITE(1,*) 'Outlier-suppression iteration',iter,' of',nter 02640 WRITE(1,*) '==========================================' 02641 WRITE(1,*) ' ' 02642 imeas=0 02643 END IF 02644 02645 ! ----- second loop ------------------------------------------------ 02646 ! accumulate normal equations for local fit and determine solution 02647 DO i=1,nalc 02648 blvec(i)=0.0D0 ! reset vector 02649 END DO 02650 DO i=1,(nalc*nalc+nalc)/2 ! GF: FIXME - not really, local parameter number... 02651 clmat(i)=0.0D0 ! (p)reset matrix 02652 END DO 02653 neq=0 02654 ndown=0 02655 nweig=0 02656 ist=isfrst(ibuf) 02657 nst=islast(ibuf) 02658 DO ! loop over measurements 02659 CALL isjajb(nst,ist,ja,jb,jsp) 02660 IF(ja == 0) EXIT 02661 rmeas=glder(ja) ! data 02662 rerr =glder(jb) ! ... and the error 02663 wght =1.0/rerr**2 ! weight from error 02664 neq=neq+1 ! count equation 02665 nweig=nweig+1 02666 resid=rmeas-localCorrections(neq) ! subtract previous fit 02667 IF(nloopn /= 1.AND.iter /= 1.AND.lhuber /= 0) THEN 02668 IF(iter <= 3) THEN 02669 IF(ABS(resid) > chuber*rerr) THEN ! down-weighting 02670 wght=wght*chuber*rerr/ABS(resid) 02671 ndown=ndown+1 02672 END IF 02673 ELSE ! Cauchy 02674 wght=wght/(1.0+(resid/rerr/cauchy)**2) 02675 END IF 02676 END IF 02677 02678 IF(lprnt.AND.iter /= 1.AND.nter /= 1) THEN 02679 chast=' ' 02680 IF(ABS(resid) > chuber*rerr) chast='* ' 02681 IF(ABS(resid) > 3.0*rerr) chast='** ' 02682 IF(ABS(resid) > 6.0*rerr) chast='***' 02683 IF(imeas == 0) WRITE(1,*) 'Second loop: accumulate' 02684 IF(imeas == 0) WRITE(1,103) 02685 imeas=imeas+1 02686 down=1.0/SQRT(wght) 02687 r1=resid/rerr 02688 r2=resid/down 02689 WRITE(1,104) imeas,rmeas,resid,rerr,r1,chast,r2 02690 END IF 02691 103 FORMAT(' index corrvalue residuum sigma', & 02692 ' nresid cnresid') 02693 104 FORMAT(i6,2X,2G12.4,' +-',g12.4,f7.2,1X,a3,f8.2) 02694 02695 DO j=1,jb-ja-1 ! normal equations, local parameter loop 02696 ij=inder(ja+j) ! local parameter index J 02697 blvec(ij)=blvec(ij)+DBLE(wght)*DBLE(rmeas)*DBLE(glder(ja+j)) 02698 DO k=1,j 02699 ik=inder(ja+k) ! local parameter index K 02700 jk=ijsym(ij,ik) ! index in symmetric matrix 02701 clmat(jk)=clmat(jk) & ! force double precision 02702 +DBLE(wght)*DBLE(glder(ja+j))*DBLE(glder(ja+k)) 02703 ! check for band matrix substructure 02704 IF (iter == 1) THEN 02705 id=IABS(ij-ik)+1 02706 im=MIN(ij,ik) 02707 ibandh(id)=MAX(ibandh(id),im) 02708 END IF 02709 END DO 02710 END DO 02711 END DO 02712 ! for non trivial fits check for bordered band matrix structure 02713 IF (iter == 1.AND.nalc > 5.AND.lfitbb > 0) THEN 02714 kx=-1 02715 kbdr=0 02716 kbdrx=0 02717 icmn=nalc**3 ! cost (*6) should improve by at least factor 2 02718 DO k=nalc,2,-1 02719 kbnd=k-2 02720 kbdr=MAX(kbdr,ibandh(k)) 02721 icost=6*(nalc-kbdr)*(kbnd+kbdr+1)**2+2*kbdr**3 02722 IF (icost < icmn) THEN 02723 icmn=icost 02724 kx=k 02725 kbdrx=kbdr 02726 END IF 02727 END DO 02728 IF (kx > 0) THEN 02729 mbnd=kx-2 02730 mbdr=kbdrx 02731 END IF 02732 END IF 02733 02734 IF (mbnd >= 0) THEN 02735 ! fast solution for border banded matrix (inverse for ICALCM>0) 02736 IF (nloopn == 1) THEN 02737 nbndr=nbndr+1 02738 nbdrx=MAX(nbdrx,mbdr) 02739 nbndx=MAX(nbndx,mbnd) 02740 END IF 02741 02742 inv=0 02743 IF (nloopn <= lfitnp.AND.iter == 1) inv=1 ! band part of inverse (for pulls) 02744 IF (icalcm == 1.OR.lprnt) inv=2 ! complete inverse 02745 CALL sqmibb(clmat,blvec,nalc,mbdr,mbnd,inv,nrank, & 02746 vbnd,vbdr,aux,vbk,vzru,scdiag,scflag) 02747 ELSE 02748 ! full inversion and solution 02749 inv=2 02750 CALL sqminv(clmat,blvec,nalc,nrank,scdiag,scflag) 02751 END IF 02752 ! check for NaNs 02753 nan=0 02754 DO k=1, nalc 02755 IF ((.NOT.(blvec(k) <= 0.0D0)).AND. (.NOT.(blvec(k) > 0.0D0))) nan=nan+1 02756 END DO 02757 02758 IF(lprnt) THEN 02759 WRITE(1,*) ' ' 02760 WRITE(1,*) 'Parameter determination:',nalc,' parameters,', ' rank=',nrank 02761 WRITE(1,*) '-----------------------' 02762 IF(ndown /= 0) WRITE(1,*) ' ',ndown,' data down-weighted' 02763 WRITE(1,*) ' ' 02764 END IF 02765 02766 ! ----- third loop ------------------------------------------------- 02767 ! calculate single residuals remaining after local fit and chi^2 02768 02769 summ=0.0D0 02770 suwt=0.0 02771 neq=0 02772 imeas=0 02773 ist=isfrst(ibuf) 02774 nst=islast(ibuf) 02775 DO ! loop over measurements 02776 CALL isjajb(nst,ist,ja,jb,jsp) 02777 IF(ja == 0) EXIT 02778 rmeas=glder(ja) ! data (global contrib. subtracted) 02779 rerr =glder(jb) ! ... and the error 02780 wght =1.0/rerr**2 ! weight from error 02781 neq=neq+1 ! count equation 02782 rmloc=0.0 ! local fit result reset 02783 DO j=1,jb-ja-1 ! local parameter loop 02784 ij=inder(ja+j) 02785 rmloc=rmloc+glder(ja+j)*SNGL(blvec(ij)) ! local fit result 02786 END DO 02787 localCorrections(neq)=rmloc ! save local fit result 02788 rmeas=rmeas-rmloc ! reduced to residual 02789 02790 ! calculate pulls? (needs covariance matrix) 02791 IF(iter == 1.AND.inv > 0.AND.nloopn <= lfitnp) THEN 02792 dvar=0.0D0 02793 DO j=1,jb-ja-1 02794 ij=inder(ja+j) 02795 DO k=1,jb-ja-1 02796 ik=inder(ja+k) 02797 jk=ijsym(ij,ik) 02798 dvar=dvar+clmat(jk)*DBLE(glder(ja+j))*DBLE(glder(ja+k)) 02799 END DO 02800 END DO 02801 ! some variance left to define a pull? 02802 IF (0.999999D0/DBLE(wght) > dvar) THEN 02803 pull=rmeas/SNGL(DSQRT(1.0D0/DBLE(wght)-dvar)) 02804 IF (lhist) THEN 02805 IF (jb < ist) THEN 02806 CALL hmpent(13,pull) ! histogram pull 02807 CALL gmpms(5,REC,pull) 02808 ELSE 02809 CALL hmpent(14,pull) ! histogram pull 02810 END IF 02811 END IF 02812 END IF 02813 END IF 02814 02815 IF(iter == 1.AND.jb < ist.AND.lhist) & 02816 CALL gmpms(4,REC,rmeas/rerr) ! residual (with global deriv.) 02817 02818 dchi2=wght*rmeas*rmeas 02819 ! DCHIT=DCHI2 02820 resid=rmeas 02821 IF(nloopn /= 1.AND.iter /= 1.AND.lhuber /= 0) THEN 02822 IF(iter <= 3) THEN 02823 IF(ABS(resid) > chuber*rerr) THEN ! down-weighting 02824 wght=wght*chuber*rerr/ABS(resid) 02825 dchi2=2.0*chuber*(ABS(resid)/rerr-0.5*chuber) 02826 END IF 02827 ELSE 02828 wght=wght/(1.0+(resid/rerr/cauchy)**2) 02829 dchi2=LOG(1.0+(resid/rerr/cauchy)**2)*cauchy**2 02830 END IF 02831 END IF 02832 02833 down=1.0/SQRT(wght) 02834 02835 ! SUWT=SUWT+DCHI2/DCHIT 02836 suwt=suwt+rerr/down 02837 IF(lprnt) THEN 02838 chast=' ' 02839 IF(ABS(resid) > chuber*rerr) chast='* ' 02840 IF(ABS(resid) > 3.0*rerr) chast='** ' 02841 IF(ABS(resid) > 6.0*rerr) chast='***' 02842 IF(imeas == 0) WRITE(1,*) 'Third loop: single residuals' 02843 IF(imeas == 0) WRITE(1,105) 02844 imeas=imeas+1 02845 r1=resid/rerr 02846 r2=resid/down 02847 IF(resid < 0.0) r1=-r1 02848 IF(resid < 0.0) r2=-r2 02849 WRITE(1,106) imeas,glder(ja),rmeas,rerr,r1,chast,r2 02850 END IF 02851 105 FORMAT(' index corrvalue residuum sigma', & 02852 ' nresid cnresid') 02853 106 FORMAT(i6,2X,2G12.4,' +-',g12.4,f7.2,1X,a3,f8.2) 02854 02855 IF(iter == nter) THEN 02856 readBufferDataF(ja)=rmeas ! store remaining residual 02857 resmax=MAX(resmax,ABS(rmeas)/rerr) 02858 END IF 02859 02860 IF(iter == 1.AND.lhist) THEN 02861 IF (jb < ist) THEN 02862 CALL hmpent( 3,rmeas/rerr) ! histogram norm residual 02863 ELSE 02864 CALL hmpent(12,rmeas/rerr) ! histogram norm residual 02865 END IF 02866 END IF 02867 summ=summ+dchi2 ! accumulate chi-square sum 02868 END DO 02869 ndf=neq-nrank 02870 ! external seed ? 02871 dch2sd=0.0D0 02872 resing=(FLOAT(nweig)-suwt)/FLOAT(nweig) 02873 IF (lhist) THEN 02874 IF(iter == 1) CALL hmpent( 5,FLOAT(ndf)) ! histogram Ndf 02875 IF(iter == 1) CALL hmpent(11,FLOAT(nalc)) ! histogram Nlocal 02876 IF(nloopn == 2.AND.iter == nter) CALL hmpent(6,resing) 02877 END IF 02878 IF(lprnt) THEN 02879 WRITE(1,*) ' ' 02880 WRITE(1,*) 'Chi^2=',summ,' at',ndf,' degrees of freedom: ', & 02881 '3-sigma limit is',chindl(3,ndf)*FLOAT(ndf) 02882 WRITE(1,*) suwt,' is sum of factors, compared to',nweig, & 02883 ' Downweight fraction:',resing 02884 END IF 02885 IF(nrank /= nalc.OR.nan > 0) THEN 02886 nrej(0)=nrej(0)+1 ! count cases 02887 IF (nrec3 == maxi4) nrec3=nrc 02888 IF(lprnt) THEN 02889 WRITE(1,*) ' rank deficit/NaN ', nalc, nrank, nan 02890 WRITE(1,*) ' ---> rejected!' 02891 END IF 02892 GO TO 90 02893 END IF 02894 IF(ndf <= 0) THEN 02895 nrej(1)=nrej(1)+1 ! count cases 02896 IF(lprnt) THEN 02897 WRITE(1,*) ' ---> rejected!' 02898 END IF 02899 GO TO 90 02900 END IF 02901 02902 chndf=SNGL(summ/dfloat(ndf)) 02903 02904 IF(iter == 1.AND.lhist) CALL hmpent(4,chndf) ! histogram chi^2/Ndf 02905 END DO ! outlier iteration loop 02906 02907 ndfs=ndfs+ndf ! (local) sum of Ndf 02908 sndf=sndf+dfloat(ndf)*dw1 ! (local) weighted sum of Ndf 02909 02910 ! ----- reject eventually ------------------------------------------ 02911 02912 IF(newite.AND.iterat == 2) THEN ! find record with largest Chi^2/Ndf 02913 IF(nrecp2 < 0.AND.chndf > writeBufferData(2,iproc+1)) THEN 02914 writeBufferData(2,iproc+1)=chndf 02915 writeBufferInfo(7,iproc+1)=nrc 02916 END IF 02917 END IF 02918 02919 chichi=chindl(3,ndf)*FLOAT(ndf) 02920 ! GF IF(SUMM.GT.50.0*CHICHI) THEN ! huge 02921 ! CHK CHICUT<0: NO cut (1st iteration) 02922 IF(chicut >= 0.0) THEN 02923 IF(summ > chhuge*chichi) THEN ! huge 02924 nrej(2)=nrej(2)+1 ! count cases with huge chi^2 02925 IF(lprnt) THEN 02926 WRITE(1,*) ' ---> rejected!' 02927 END IF 02928 GO TO 90 02929 END IF 02930 02931 IF(chicut > 0.0) THEN 02932 chlimt=chicut*chichi 02933 ! WRITE(*,*) 'chi^2 ',SUMM,CHLIMT,CHICUT,CHINDL(3,NDF),NDF 02934 IF(summ > chlimt) THEN 02935 IF(lprnt) THEN 02936 WRITE(1,*) ' ---> rejected!' 02937 END IF 02938 ! add to FVALUE 02939 dchi2=chlimt ! total contribution limit 02940 dchi2s=dchi2s+dchi2*dw1 ! add total contribution 02941 nrej(3)=nrej(3)+1 ! count cases with large chi^2 02942 GO TO 90 02943 END IF 02944 END IF 02945 END IF 02946 02947 IF(lhuber > 1.AND.dwcut /= 0.0.AND.resing > dwcut) THEN 02948 ! add to FVALUE 02949 dchi2=summ ! total contribution 02950 dchi2s=dchi2s+dchi2*dw1 ! add total contribution 02951 nrej(3)=nrej(3)+1 ! count cases with large chi^2 02952 ! WRITE(*,*) 'Downweight fraction cut ',RESING,DWCUT,SUMM 02953 IF(lprnt) THEN 02954 WRITE(1,*) ' ---> rejected!' 02955 END IF 02956 GO TO 90 02957 END IF 02958 02959 IF(newite.AND.iterat == 2) THEN ! find record with largest residual 02960 IF(nrecpr < 0.AND.resmax > writeBufferData(1,iproc+1)) THEN 02961 writeBufferData(1,iproc+1)=resmax 02962 writeBufferInfo(6,iproc+1)=nrc 02963 END IF 02964 END IF 02965 ! 'track quality' per binary file: accepted records 02966 naccf(kfl)=naccf(kfl)+1 02967 ndff(kfl) =ndff(kfl) +ndf 02968 chi2f(kfl)=chi2f(kfl)+chndf 02969 02970 ! ----- fourth loop ------------------------------------------------ 02971 ! update of global matrix and vector according to the "Millepede" 02972 ! principle, from the global/local information 02973 02974 dchi2s=dchi2s+dch2sd 02975 02976 ist=isfrst(ibuf) 02977 nst=islast(ibuf) 02978 DO ! loop over measurements 02979 CALL isjajb(nst,ist,ja,jb,jsp) 02980 IF(ja <= 0) EXIT 02981 02982 rmeas=glder(ja) ! data residual 02983 rerr =glder(jb) ! ... and the error 02984 wght =1.0/rerr**2 ! weight from measurement error 02985 dchi2=wght*rmeas*rmeas ! least-square contribution 02986 02987 IF(nloopn /= 1.AND.lhuber /= 0) THEN ! check residual 02988 resid=ABS(rmeas) 02989 IF(resid > chuber*rerr) THEN 02990 wght=wght*chuber*rerr/resid ! down-weighting 02991 dchi2=2.0*chuber*(resid/rerr-0.5*chuber) ! modified contribution 02992 END IF 02993 END IF 02994 dchi2s=dchi2s+dchi2*dw1 ! add to total objective function 02995 02996 ! global-global matrix contribution: add directly to gg-matrix 02997 02998 DO j=1,ist-jb 02999 ivgbj=globalParLabelIndex(2,inder(jb+j)) ! variable-parameter index 03000 IF(ivgbj > 0) THEN 03001 globalVector(ioffb+ivgbj)=globalVector(ioffb+ivgbj) & 03002 +dw1*wght*rmeas*glder(jb+j) ! vector !!! reverse 03003 IF(icalcm == 1) THEN 03004 ije=backIndexUsage(ioffe+ivgbj) ! get index of index, non-zero 03005 DO k=1,j 03006 ivgbk=globalParLabelIndex(2,inder(jb+k)) 03007 IF(ivgbk > 0) THEN 03008 ike=backIndexUsage(ioffe+ivgbk) ! get index of index, non-zero 03009 ia=MAX(ije,ike) ! larger 03010 ib=MIN(ije,ike) ! smaller 03011 ij=ib+(ia*ia-ia)/2 03012 writeBufferUpdates(ioffd+ij)=writeBufferUpdates(ioffd+ij) & 03013 -dw1*wght*glder(jb+j)*glder(jb+k) 03014 END IF 03015 END DO 03016 END IF 03017 END IF 03018 END DO 03019 03020 ! normal equations - rectangular matrix for global/local pars 03021 ! global-local matrix contribution: accumulate rectangular matrix 03022 IF (icalcm /= 1) CYCLE 03023 DO j=1,ist-jb 03024 ivgbj=globalParLabelIndex(2,inder(jb+j)) ! variable-parameter index 03025 IF(ivgbj > 0) THEN 03026 ije=backIndexUsage(ioffe+ivgbj) ! get index of index, non-zero 03027 DO k=1,jb-ja-1 03028 ik=inder(ja+k) ! local index 03029 jk=ik+(ije-1)*nalc ! matrix index 03030 localGlobalMatrix(jk)=localGlobalMatrix(jk)+dw2*wght*glder(jb+j)*glder(ja+k) 03031 END DO 03032 END IF 03033 END DO 03034 END DO 03035 03036 03037 ! ----- final matrix update ---------------------------------------- 03038 ! update global matrices and vectors 03039 IF(icalcm /= 1) GO TO 90 ! matrix update 03040 ! (inverse local matrix) * (rectang. matrix) -> CORM 03041 ! T 03042 ! resulting symmetrix matrix = G * Gamma^{-1} * G 03043 03044 CALL dbavat(clmat,localGlobalMatrix,writeBufferUpdates(ioffd+1),nalc,-nalg) 03045 03046 ! (rectang. matrix) * (local param vector) -> CORV 03047 ! resulting vector = G * q (q = local parameter) 03048 ! CALL DBGAX(DQ(IGLMA/2+1),BLVEC,DQ(ICORV/2+1),NALG,NALC) ! not done 03049 ! the vector update is not done, because after local fit it is zero! 03050 03051 ! update cache status 03052 writeBufferInfo(1,iproc+1)=writeBufferInfo(1,iproc+1)+1 03053 writeBufferInfo(2,iproc+1)=writeBufferInfo(2,iproc+1)+ngg 03054 writeBufferInfo(3,iproc+1)=writeBufferInfo(3,iproc+1)+nalg+2 03055 ! check free space 03056 nfred=writeBufferHeader(-1)-writeBufferInfo(2,iproc+1)-writeBufferHeader(-2) 03057 nfrei=writeBufferHeader(1)-writeBufferInfo(3,iproc+1)-writeBufferHeader(2) 03058 IF (nfred < 0.OR.nfrei < 0) THEN ! need to flush 03059 nb=writeBufferInfo(1,iproc+1) 03060 joffd=writeBufferHeader(-1)*iproc ! offset data 03061 joffi=writeBufferHeader(1)*iproc+2 ! offset indices 03062 used=FLOAT(writeBufferInfo(2,iproc+1))/FLOAT(writeBufferHeader(-1)) 03063 writeBufferInfo(4,iproc+1)=writeBufferInfo(4,iproc+1) +IFIX(1000.0*used+0.5) 03064 used=FLOAT(writeBufferInfo(3,iproc+1))/FLOAT(writeBufferHeader(1)) 03065 writeBufferInfo(5,iproc+1)=writeBufferInfo(5,iproc+1) +IFIX(1000.0*used+0.5) 03066 !$OMP CRITICAL 03067 writeBufferHeader(-4)=writeBufferHeader(-4)+1 03068 writeBufferHeader(4)=writeBufferHeader(4)+1 03069 03070 DO ib=1,nb 03071 ijn=0 03072 DO in=1,writeBufferIndices(joffi) 03073 i=writeBufferIndices(joffi+in) 03074 ! DQ(IGVEC/2+I)=DQ(IGVEC/2+I)+DQ(ICORV/2+IN) ! not done: = zero 03075 DO jn=1,in 03076 ijn=ijn+1 03077 j=writeBufferIndices(joffi+jn) 03078 CALL mupdat(i,j,-writeBufferUpdates(joffd+ijn)) ! matrix update 03079 END DO 03080 END DO 03081 joffd=joffd+ijn 03082 joffi=joffi+writeBufferIndices(joffi)+2 03083 END DO 03084 !$OMP END CRITICAL 03085 ! reset counter, pointers 03086 DO k=1,3 03087 writeBufferInfo(k,iproc+1)=0 03088 END DO 03089 END IF 03090 03091 90 IF(lprnt) THEN 03092 WRITE(1,*) ' ' 03093 WRITE(1,*) '------------------ End of printout for record',nrc 03094 WRITE(1,*) ' ' 03095 END IF 03096 03097 DO i=1,nalg ! reset global index array 03098 iext=globalIndexUsage(ioffc+i) 03099 backIndexUsage(ioffe+iext)=0 03100 END DO 03101 03102 END DO 03103 !$OMP END PARALLEL DO 03104 03105 IF (icalcm == 1) THEN 03106 ! flush remaining matrices 03107 DO k=1,mthrd ! update statistics 03108 writeBufferHeader(-3)=writeBufferHeader(-3)+1 03109 used=FLOAT(writeBufferInfo(2,k))/FLOAT(writeBufferHeader(-1)) 03110 writeBufferInfo(4,k)=writeBufferInfo(4,k)+IFIX(1000.0*used+0.5) 03111 writeBufferHeader(-5)=writeBufferHeader(-5)+writeBufferInfo(4,k) 03112 writeBufferHeader(-6)=MAX(writeBufferHeader(-6),writeBufferInfo(4,k)) 03113 writeBufferInfo(4,k)=0 03114 writeBufferHeader(3)=writeBufferHeader(3)+1 03115 used=FLOAT(writeBufferInfo(3,k))/FLOAT(writeBufferHeader(1)) 03116 writeBufferInfo(5,k)=writeBufferInfo(5,k)+IFIX(1000.0*used+0.5) 03117 writeBufferHeader(5)=writeBufferHeader(5)+writeBufferInfo(5,k) 03118 writeBufferHeader(6)=MAX(writeBufferHeader(6),writeBufferInfo(5,k)) 03119 writeBufferInfo(5,k)=0 03120 END DO 03121 03122 !$OMP PARALLEL & 03123 !$OMP DEFAULT(PRIVATE) & 03124 !$OMP SHARED(writeBufferHeader,writeBufferInfo,writeBufferIndices,writeBufferUpdates,MTHRD) 03125 iproc=0 03126 !$ IPROC=OMP_GET_THREAD_NUM() ! thread number 03127 DO jproc=0,mthrd-1 03128 nb=writeBufferInfo(1,jproc+1) 03129 ! print *, ' flush end ', JPROC, NRC, NB 03130 joffd=writeBufferHeader(-1)*jproc ! offset data 03131 joffi=writeBufferHeader(1)*jproc+2 ! offset indices 03132 DO ib=1,nb 03133 ! print *, ' buf end ', JPROC,IB,writeBufferIndices(JOFFI-1),writeBufferIndices(JOFFI) 03134 ijn=0 03135 DO in=1,writeBufferIndices(joffi) 03136 i=writeBufferIndices(joffi+in) 03137 !$ IF (MOD(I,MTHRD).EQ.IPROC) THEN 03138 DO jn=1,in 03139 ijn=ijn+1 03140 j=writeBufferIndices(joffi+jn) 03141 CALL mupdat(i,j,-writeBufferUpdates(joffd+ijn)) ! matrix update 03142 END DO 03143 !$ ELSE 03144 !$ IJN=IJN+IN 03145 !$ ENDIF 03146 END DO 03147 joffd=joffd+ijn 03148 joffi=joffi+writeBufferIndices(joffi)+2 03149 END DO 03150 END DO 03151 !$OMP END PARALLEL 03152 END IF 03153 03154 IF(newite.AND.iterat == 2) THEN ! get worst records (for printrecord -1 -1) 03155 IF (nrecpr < 0) THEN 03156 DO k=1,mthrd 03157 IF (writeBufferData(1,k) > value1) THEN 03158 value1=writeBufferData(1,k) 03159 nrec1 =writeBufferInfo(6,k) 03160 END IF 03161 END DO 03162 END IF 03163 IF (nrecp2 < 0) THEN 03164 DO k=1,mthrd 03165 IF (writeBufferData(2,k) > value2) THEN 03166 value2=writeBufferData(2,k) 03167 nrec2 =writeBufferInfo(7,k) 03168 END IF 03169 END DO 03170 END IF 03171 END IF 03172 03173 END SUBROUTINE loopbf 03174 03175 03176 03177 03178 !*********************************************************************** 03179 03191 SUBROUTINE prtglo 03192 USE mpmod 03193 03194 IMPLICIT NONE 03195 REAL:: dpa 03196 REAL:: err 03197 REAL:: gcor 03198 INTEGER:: i 03199 INTEGER:: ie 03200 INTEGER:: iev 03201 INTEGER:: ij 03202 INTEGER:: imin 03203 INTEGER:: iprlim 03204 INTEGER:: isub 03205 INTEGER:: itgbi 03206 INTEGER:: itgbl 03207 INTEGER:: ivgbi 03208 INTEGER:: j 03209 INTEGER:: label 03210 INTEGER:: lup 03211 REAL:: par 03212 03213 DOUBLE PRECISION:: diag 03214 DOUBLE PRECISION::gmati 03215 DOUBLE PRECISION::gcor2 03216 INTEGER:: labele(3) 03217 INTEGER(KIND=large):: ii 03218 REAL:: compnt(3) 03219 SAVE 03220 ! ... 03221 03222 lup=09 03223 CALL mvopen(lup,'millepede.res') 03224 03225 WRITE(*,*) ' ' 03226 WRITE(*,*) ' Result of fit for global parameters' 03227 WRITE(*,*) ' ===================================' 03228 WRITE(*,*) ' ' 03229 03230 WRITE(*,101) 03231 03232 WRITE(lup,*) 'Parameter ! first 3 elements per line are', & 03233 ' significant (if used as input)' 03234 iprlim=10 03235 DO itgbi=1,ntgb ! all parameter variables 03236 itgbl=globalParLabelIndex(1,itgbi) 03237 ivgbi=globalParLabelIndex(2,itgbi) 03238 par=SNGL(globalParameter(itgbi)) ! initial value 03239 IF(ivgbi > 0) THEN 03240 dpa=par-globalParStart(itgbi) ! difference 03241 IF(metsol == 1.OR.metsol == 2) THEN 03242 ii=ivgbi 03243 ii=(ii*ii+ii)/2 03244 gmati=globalMatD(ii) 03245 ERR=SQRT(ABS(SNGL(gmati))) 03246 IF(gmati < 0.0D0) ERR=-ERR 03247 diag=workspaceD(ivgbi) 03248 gcor=-1.0 03249 IF(gmati*diag > 0.0D0) THEN ! global correlation 03250 gcor2=1.0D0-1.0D0/(gmati*diag) 03251 IF(gcor2 >= 0.0D0.AND.gcor2 <= 1.0D0) gcor=SNGL(DSQRT(gcor2)) 03252 END IF 03253 END IF 03254 END IF 03255 IF(itgbi <= iprlim) THEN 03256 IF(ivgbi <= 0) THEN 03257 WRITE(* ,102) itgbl,par,globalParPreSigma(itgbi) 03258 ELSE 03259 IF(metsol == 1.OR.metsol == 2) THEN 03260 IF (igcorr == 0) THEN 03261 WRITE(*,102) itgbl,par,globalParPreSigma(itgbi),dpa,ERR 03262 ELSE 03263 WRITE(*,102) itgbl,par,globalParPreSigma(itgbi),dpa,ERR,gcor 03264 END IF 03265 ELSE 03266 WRITE(*,102) itgbl,par,globalParPreSigma(itgbi),dpa 03267 END IF 03268 END IF 03269 ELSE IF(itgbi == iprlim+1) THEN 03270 WRITE(* ,*) '... (further printout suppressed, but see log file)' 03271 END IF 03272 03273 ! file output 03274 IF(ivgbi <= 0) THEN 03275 WRITE(lup,102) itgbl,par,globalParPreSigma(itgbi) 03276 ELSE 03277 IF(metsol == 1.OR.metsol == 2) THEN 03278 IF (igcorr == 0) THEN 03279 WRITE(lup,102) itgbl,par,globalParPreSigma(itgbi),dpa,ERR 03280 ELSE 03281 WRITE(lup,102) itgbl,par,globalParPreSigma(itgbi),dpa,ERR,gcor 03282 END IF 03283 ELSE 03284 WRITE(lup,102) itgbl,par,globalParPreSigma(itgbi),dpa 03285 END IF 03286 END IF 03287 END DO 03288 REWIND lup 03289 CLOSE(UNIT=lup) 03290 03291 IF(metsol == 2) THEN ! diagonalisation: write eigenvectors 03292 CALL mvopen(lup,'millepede.eve') 03293 imin=1 03294 DO i=nagb,1,-1 03295 IF(workspaceEigenValues(i) > 0.0D0) THEN 03296 imin=i ! index of smallest pos. eigenvalue 03297 EXIT 03298 ENDIF 03299 END DO 03300 iev=0 03301 03302 DO isub=0,MIN(15,imin-1) 03303 IF(isub < 10) THEN 03304 i=imin-isub 03305 ELSE 03306 i=isub-9 03307 END IF 03308 03309 ! DO I=IMIN,MAX(1,IMIN-9),-1 ! backward loop, up to 10 vectors 03310 WRITE(*,*) 'Eigenvector ',i,' with eigenvalue',workspaceEigenValues(i) 03311 WRITE(lup,*) 'Eigenvector ',i,' with eigenvalue',workspaceEigenValues(i) 03312 DO j=1,nagb 03313 ij=j+(i-1)*nagb ! index with eigenvector array 03314 IF(j <= nvgb) THEN 03315 itgbi=globalParVarToTotal(j) 03316 label=globalParLabelIndex(1,itgbi) 03317 ELSE 03318 label=nvgb-j ! label negative for constraints 03319 END IF 03320 iev=iev+1 03321 labele(iev)=label 03322 compnt(iev)=SNGL(workspaceEigenVectors(ij)) ! component 03323 IF(iev == 3) THEN 03324 WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev) 03325 iev=0 03326 END IF 03327 END DO 03328 IF(iev /= 0) WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev) 03329 iev=0 03330 WRITE(lup,*) ' ' 03331 END DO 03332 03333 END IF 03334 03335 101 FORMAT(1X,' label parameter presigma differ', & 03336 ' error'/ 1X,'-----------',4X,4('-------------')) 03337 102 FORMAT(i10,2X,4G14.5,f8.3) 03338 103 FORMAT(3(i11,f11.7,2X)) 03339 END SUBROUTINE prtglo ! print final log file 03340 03341 03350 03351 SUBROUTINE avprod(n,x,b) 03352 USE mpmod 03353 03354 IMPLICIT NONE 03355 INTEGER :: i 03356 INTEGER :: iencdb 03357 INTEGER :: iencdm 03358 INTEGER :: iproc 03359 INTEGER :: ir 03360 INTEGER :: j 03361 INTEGER :: jc 03362 INTEGER :: jj 03363 INTEGER :: jn 03364 03365 INTEGER, INTENT(IN) :: n 03366 DOUBLE PRECISION, INTENT(IN) :: x(n) 03367 DOUBLE PRECISION, INTENT(OUT) :: b(n) 03368 INTEGER(KIND=large) :: k 03369 INTEGER(KIND=large) :: kk 03370 INTEGER(KIND=large) :: kl 03371 INTEGER(KIND=large) :: ku 03372 INTEGER(KIND=large) :: ll 03373 INTEGER(KIND=large) :: lj 03374 INTEGER(KIND=large) :: indij 03375 INTEGER(KIND=large) :: indid 03376 INTEGER(KIND=large) :: ij 03377 INTEGER :: ichunk 03378 !$ INTEGER OMP_GET_THREAD_NUM 03379 SAVE 03380 ! ... 03381 !$ DO i=1,n 03382 !$ b(i)=0.0D0 ! reset 'global' B() 03383 !$ END DO 03384 ichunk=MIN((n+mthrd-1)/mthrd/8+1,1024) 03385 IF(matsto == 1) THEN 03386 ! full symmetric matrix 03387 ! parallelize row loop 03388 ! private copy of B(N) for each thread, combined at end, init with 0. 03389 ! slot of 1024 'I' for next idle thread 03390 !$OMP PARALLEL DO & 03391 !$OMP PRIVATE(J,IJ) & 03392 !$OMP REDUCTION(+:B) & 03393 !$OMP SCHEDULE(DYNAMIC,ichunk) 03394 DO i=1,n 03395 ij=i 03396 ij=(ij*ij-ij)/2 03397 b(i)=globalMatD(ij+i)*x(i) 03398 DO j=1,i-1 03399 b(j)=b(j)+globalMatD(ij+j)*x(i) 03400 b(i)=b(i)+globalMatD(ij+j)*x(j) 03401 END DO 03402 END DO 03403 !$OMP END PARALLEL DO 03404 ELSE 03405 ! sparse, compressed matrix 03406 IF(sparseMatrixOffsets(2,1) /= n+1) STOP 'AVPROD: mismatched vector and matrix' 03407 iencdb=nencdb 03408 iencdm=ishft(1,iencdb)-1 03409 ! parallelize row loop 03410 ! slot of 1024 'I' for next idle thread 03411 !$OMP PARALLEL DO & 03412 !$OMP PRIVATE(IR,K,KK,LL,KL,KU,INDID,INDIJ,J,JC,JN,LJ,JJ) & 03413 !$OMP REDUCTION(+:B) & 03414 !$OMP SCHEDULE(DYNAMIC,ichunk) 03415 DO i=1,n 03416 iproc=0 03417 !$ IPROC=OMP_GET_THREAD_NUM() ! thread number 03418 b(i)=globalMatD(i)*x(i) ! diagonal elements 03419 ! ! off-diagonals double precision 03420 ir=i 03421 kk=sparseMatrixOffsets(1,ir) ! offset in 'd' (column lists) 03422 ll=sparseMatrixOffsets(2,ir) ! offset in 'j' (matrix) 03423 kl=0 03424 ku=sparseMatrixOffsets(1,ir+1)-1-kk 03425 indid=kk 03426 indij=ll 03427 IF (sparseMatrixColumns(indid) /= 0) THEN ! no compression 03428 DO k=kl,ku 03429 j=sparseMatrixColumns(indid+k) 03430 b(j)=b(j)+globalMatD(indij+k)*x(i) 03431 b(i)=b(i)+globalMatD(indij+k)*x(j) 03432 END DO 03433 ELSE 03434 lj=0 03435 ku=((ku+1)*8)/9-1 ! number of regions (-1) 03436 indid=indid+ku/8+1 ! skip group offsets 03437 DO kl=0,ku 03438 jc=sparseMatrixColumns(indid+kl) 03439 j=ishft(jc,-iencdb) 03440 jn=IAND(jc, iencdm) 03441 DO jj=1,jn 03442 b(j)=b(j)+globalMatD(indij+lj)*x(i) 03443 b(i)=b(i)+globalMatD(indij+lj)*x(j) 03444 j=j+1 03445 lj=lj+1 03446 END DO 03447 END DO 03448 END IF 03449 03450 IF (nspc > 1) THEN 03451 ir=i+n+1 ! off-diagonals single precision 03452 kk=sparseMatrixOffsets(1,ir) ! offset in 'd' (column lists) 03453 ll=sparseMatrixOffsets(2,ir) ! offset in '.' (matrix) 03454 kl=0 03455 ku=sparseMatrixOffsets(1,ir+1)-1-kk 03456 indid=kk 03457 indij=ll 03458 IF (sparseMatrixColumns(indid) /= 0) THEN ! no compression 03459 DO k=kl,ku 03460 j=sparseMatrixColumns(indid+k) 03461 b(j)=b(j)+DBLE(globalMatF(indij+k))*x(i) 03462 b(i)=b(i)+DBLE(globalMatF(indij+k))*x(j) 03463 END DO 03464 ELSE 03465 lj=0 03466 ku=((ku+1)*8)/9-1 ! number of regions (-1) 03467 indid=indid+ku/8+1 ! skip group offsets 03468 DO kl=0,ku 03469 jc=sparseMatrixColumns(indid+kl) 03470 j=ishft(jc,-iencdb) 03471 jn=IAND(jc, iencdm) 03472 DO jj=1,jn 03473 b(j)=b(j)+DBLE(globalMatF(indij+lj))*x(i) 03474 b(i)=b(i)+DBLE(globalMatF(indij+lj))*x(j) 03475 j=j+1 03476 lj=lj+1 03477 END DO 03478 END DO 03479 END IF 03480 END IF 03481 END DO 03482 !$OMP END PARALLEL DO 03483 ENDIF 03484 03485 END SUBROUTINE avprod 03486 03494 03495 FUNCTION ijadd(itema,itemb) ! index using "d" and "z" 03496 USE mpmod 03497 03498 IMPLICIT NONE 03499 INTEGER :: iencdb 03500 INTEGER :: iencdm 03501 INTEGER :: isgn 03502 INTEGER :: ispc 03503 INTEGER :: item2 03504 INTEGER :: jtem 03505 INTEGER :: jtemc 03506 INTEGER :: jtemn 03507 03508 INTEGER, INTENT(IN) :: itema 03509 INTEGER, INTENT(IN) :: itemb 03510 03511 INTEGER(KIND=large) :: ijadd 03512 INTEGER(KIND=large) :: k 03513 INTEGER(KIND=large) :: kk 03514 INTEGER(KIND=large) :: kl 03515 INTEGER(KIND=large) :: ku 03516 INTEGER(KIND=large) :: indid 03517 INTEGER(KIND=large) :: nd 03518 INTEGER(KIND=large) :: ll 03519 INTEGER(KIND=large) :: k8 03520 INTEGER(KIND=large) :: item1 03521 ! ... 03522 ijadd=0 03523 nd=sparseMatrixOffsets(2,1)-1 ! dimension of matrix 03524 item1=MAX(itema,itemb) ! larger index 03525 item2=MIN(itema,itemb) ! smaller index 03526 IF(item2 <= 0.OR.item1 > nd) RETURN 03527 IF(item1 == item2) THEN ! diagonal element 03528 ijadd=item1 03529 RETURN 03530 END IF 03531 ! ! off-diagonal element 03532 iencdb=nencdb ! encoding info 03533 iencdm=ishft(1,iencdb)-1 03534 isgn=-1 03535 outer: DO ispc=1,nspc 03536 kk=sparseMatrixOffsets(1,item1) ! offset in 'd' (column lists) 03537 ll=sparseMatrixOffsets(2,item1) ! offset in 'j' (matrix) 03538 kl=0 03539 ku=sparseMatrixOffsets(1,item1+1)-1-kk 03540 indid=kk 03541 item1=item1+nd+1 03542 isgn=-isgn 03543 IF (sparseMatrixColumns(indid) == 0) THEN ! compression ? 03544 03545 ku=((ku+1)*8)/9-1 ! number of regions (-1) 03546 indid=indid+ku/8+1 ! skip group offsets 03547 kl=0 03548 IF(ku < kl) CYCLE outer ! not found 03549 DO 03550 k=(kl+ku)/2 ! binary search 03551 jtemc=sparseMatrixColumns(indid+k) ! compressed information 03552 jtem =ishft(jtemc,-iencdb) ! first column of region 03553 jtemn=jtem+IAND(jtemc,iencdm) ! first column after region 03554 IF(item2 >= jtem.AND.item2 < jtemn) EXIT ! found 03555 IF(item2 < jtem) THEN 03556 ku=k-1 03557 ELSE IF(item2 >= jtemn) THEN 03558 kl=k+1 03559 END IF 03560 IF(kl <= ku) CYCLE 03561 CYCLE outer ! not found 03562 END DO 03563 k8=k/8 ! region group (-1) 03564 ll=ll+sparseMatrixColumns(kk+k8) ! offset for group of (8) regions 03565 DO kl=k8*8,k-1 03566 ll=ll+IAND(sparseMatrixColumns(indid+kl),iencdm) ! add region lengths 03567 END DO 03568 ijadd=ll+item2-jtem 03569 03570 ELSE 03571 03572 IF(ku < kl) CYCLE outer ! not found 03573 DO 03574 k=(kl+ku)/2 ! binary search 03575 jtem=sparseMatrixColumns(indid+k) 03576 jtemn=jtem 03577 IF(item2 == jtem) EXIT ! found 03578 IF(item2 < jtem) THEN 03579 ku=k-1 03580 ELSE IF(item2 > jtem) THEN 03581 kl=k+1 03582 END IF 03583 IF(kl <= ku) CYCLE 03584 CYCLE outer ! not found 03585 END DO 03586 ijadd=ll+k 03587 03588 END IF 03589 ijadd=ijadd*isgn 03590 RETURN 03591 END DO outer 03592 03593 END FUNCTION ijadd 03594 03603 03604 SUBROUTINE sechms(deltat,nhour,minut,secnd) 03605 IMPLICIT NONE 03606 REAL, INTENT(IN) :: deltat 03607 INTEGER, INTENT(OUT) :: minut 03608 INTEGER, INTENT(OUT):: nhour 03609 REAL, INTENT(OUT):: secnd 03610 INTEGER:: nsecd 03611 ! DELTAT = time in sec -> NHOUR,MINUT,SECND 03612 ! ... 03613 nsecd=nint(deltat) ! -> integer 03614 nhour=nsecd/3600 03615 minut=nsecd/60-60*nhour 03616 secnd=deltat-60*(minut+60*nhour) 03617 END SUBROUTINE sechms 03618 03646 03647 INTEGER FUNCTION inone(item) ! translate 1-D identifier to nrs 03648 USE mpmod 03649 USE mpdalc 03650 03651 IMPLICIT NONE 03652 INTEGER, INTENT(IN) :: item 03653 INTEGER :: j 03654 INTEGER :: k 03655 INTEGER :: iprime 03656 INTEGER(KIND=large) :: length 03657 INTEGER(KIND=large), PARAMETER :: two = 2 03658 03659 inone=0 03660 IF(item <= 0) RETURN 03661 IF(globalParHeader(-1) == 0) THEN 03662 length=128 ! initial number 03663 CALL mpalloc(globalParLabelIndex,two,length,'INONE: label & index') 03664 CALL mpalloc(globalParHashTable,2*length,'INONE: hash pointer') 03665 globalParHashTable = 0 03666 globalParHeader(-0)=INT(length) ! length of labels/indices 03667 globalParHeader(-1)=0 ! number of stored items 03668 globalParHeader(-2)=0 ! =0 during build-up 03669 globalParHeader(-3)=INT(length) ! next number 03670 globalParHeader(-4)=iprime(globalParHeader(-0)) ! prime number 03671 globalParHeader(-5)=0 ! number of overflows 03672 globalParHeader(-6)=0 ! nr of variable parameters 03673 END IF 03674 outer: DO 03675 j=1+MOD(item,globalParHeader(-4))+globalParHeader(-0) 03676 inner: DO ! normal case: find item 03677 k=j 03678 j=globalParHashTable(k) 03679 IF(j == 0) EXIT inner ! unused hash code 03680 IF(item == globalParLabelIndex(1,j)) EXIT outer ! found 03681 END DO inner 03682 ! not found 03683 IF(globalParHeader(-1) == globalParHeader(-0).OR.globalParHeader(-2) /= 0) THEN 03684 globalParHeader(-5)=globalParHeader(-5)+1 ! overflow 03685 j=0 03686 RETURN 03687 END IF 03688 globalParHeader(-1)=globalParHeader(-1)+1 ! increase number of elements 03689 globalParHeader(-3)=globalParHeader(-1) 03690 j=globalParHeader(-1) 03691 globalParHashTable(k)=j ! hash index 03692 globalParLabelIndex(1,j)=item ! add new item 03693 globalParLabelIndex(2,j)=0 ! reset counter 03694 IF(globalParHeader(-1) /= globalParHeader(-0)) EXIT outer 03695 ! update with larger dimension and redefine index 03696 globalParHeader(-3)=globalParHeader(-3)*2 03697 CALL upone 03698 IF (lvllog > 1) WRITE(lunlog,*) 'INONE: array increased to', & 03699 globalParHeader(-3),' words' 03700 END DO outer 03701 03702 IF(globalParHeader(-2) == 0) THEN 03703 globalParLabelIndex(2,j)=globalParLabelIndex(2,j)+1 ! increase counter 03704 globalParHeader(-7)=globalParHeader(-7)+1 03705 END IF 03706 inone=j 03707 END FUNCTION inone 03708 03710 SUBROUTINE upone 03711 USE mpmod 03712 USE mpdalc 03713 03714 IMPLICIT NONE 03715 INTEGER :: i 03716 INTEGER :: j 03717 INTEGER :: k 03718 INTEGER :: iprime 03719 INTEGER :: nused 03720 LOGICAL :: finalUpdate 03721 INTEGER(KIND=large) :: oldLength 03722 INTEGER(KIND=large) :: newLength 03723 INTEGER(KIND=large), PARAMETER :: two = 2 03724 INTEGER, DIMENSION(:,:), ALLOCATABLE :: tempArr 03725 SAVE 03726 ! ... 03727 finalUpdate=(globalParHeader(-3) == globalParHeader(-1)) 03728 IF(finalUpdate) THEN ! final (cleanup) call 03729 CALL sort2k(globalParLabelIndex,globalParHeader(-1)) ! sort items 03730 END IF 03731 ! save old LabelIndex 03732 nused = globalParHeader(-1) 03733 oldLength = globalParHeader(-0) 03734 CALL mpalloc(tempArr,two,oldLength,'INONE: temp array') 03735 tempArr(:,1:nused)=globalParLabelIndex(:,1:nused) 03736 CALL mpdealloc(globalParLabelIndex) 03737 CALL mpdealloc(globalParHashTable) 03738 ! create new LabelIndex 03739 newLength = globalParHeader(-3) 03740 CALL mpalloc(globalParLabelIndex,two,newLength,'INONE: label & index') 03741 CALL mpalloc(globalParHashTable,2*newLength,'INONE: hash pointer') 03742 globalParHashTable = 0 03743 globalParLabelIndex(:,1:nused) = tempArr(:,1:nused) ! copy back saved content 03744 CALL mpdealloc(tempArr) 03745 globalParHeader(-0)=INT(newLength) ! length of labels/indices 03746 globalParHeader(-3)=globalParHeader(-1) 03747 globalParHeader(-4)=iprime(globalParHeader(-0)) ! prime number < LNDA 03748 ! redefine hash 03749 outer: DO i=1,globalParHeader(-1) 03750 j=1+MOD(globalParLabelIndex(1,i),globalParHeader(-4))+globalParHeader(-0) 03751 inner: DO 03752 k=j 03753 j=globalParHashTable(k) 03754 IF(j == 0) EXIT inner ! unused hash code 03755 IF(j == i) CYCLE outer ! found 03756 ENDDO inner 03757 globalParHashTable(k)=i 03758 END DO outer 03759 IF(.NOT.finalUpdate) RETURN 03760 03761 globalParHeader(-2)=1 ! set flag to inhibit further updates 03762 IF (lvllog > 1) THEN 03763 WRITE(lunlog,*) ' ' 03764 WRITE(lunlog,*) 'INONE: array reduced to',newLength,' words' 03765 WRITE(lunlog,*) 'INONE:',globalParHeader(-1),' items stored.' 03766 END IF 03767 END SUBROUTINE upone ! update, redefine 03768 03773 03774 INTEGER FUNCTION iprime(n) 03775 IMPLICIT NONE 03776 INTEGER, INTENT(IN) :: n 03777 INTEGER :: nprime 03778 INTEGER :: nsqrt 03779 INTEGER :: i 03780 ! ... 03781 SAVE 03782 nprime=n ! max number 03783 IF(MOD(nprime,2) == 0) nprime=nprime+1 ! ... odd number 03784 outer: DO 03785 nprime=nprime-2 ! next lower odd number 03786 nsqrt=IFIX(SQRT(FLOAT(nprime))) 03787 DO i=3,nsqrt,2 ! 03788 IF(i*(nprime/i) == nprime) CYCLE outer ! test prime number 03789 END DO 03790 EXIT outer ! found 03791 END DO outer 03792 iprime=nprime 03793 END FUNCTION iprime 03794 03804 SUBROUTINE loop1 03805 USE mpmod 03806 USE mpdalc 03807 03808 IMPLICIT NONE 03809 INTEGER :: i 03810 INTEGER :: idum 03811 INTEGER :: in 03812 INTEGER :: indab 03813 INTEGER :: itgbi 03814 INTEGER :: itgbl 03815 INTEGER :: ivgbi 03816 INTEGER :: j 03817 INTEGER :: mqi 03818 INTEGER :: nc21 03819 INTEGER :: nr 03820 INTEGER :: nwrd 03821 INTEGER :: inone 03822 REAL :: param 03823 REAL :: presg 03824 REAL :: prewt 03825 03826 REAL :: plvs(3) ! vector array: real and ... 03827 INTEGER :: lpvs(3) ! ... integer 03828 EQUIVALENCE (plvs(1),lpvs(1)) 03829 INTEGER (KIND=large) :: length 03830 SAVE 03831 ! ... 03832 WRITE(lunlog,*) ' ' 03833 WRITE(lunlog,*) 'LOOP1: starting' 03834 CALL mstart('LOOP1') 03835 ! add labels from parameter, constraints, measurements ------------- 03836 DO i=1, lenParameters 03837 idum=inone(listParameters(i)%label) 03838 END DO 03839 DO i=1, lenPreSigmas 03840 idum=inone(listPreSigmas(i)%label) 03841 END DO 03842 DO i=1, lenConstraints 03843 idum=inone(listConstraints(i)%label) 03844 END DO 03845 DO i=1, lenMeasurements 03846 idum=inone(listMeasurements(i)%label) 03847 END DO 03848 03849 IF(globalParHeader(-1) /= 0) THEN 03850 WRITE(lunlog,*) 'LOOP1:',globalParHeader(-1), ' labels from txt data stored' 03851 END IF 03852 WRITE(lunlog,*) 'LOOP1: reading data files' 03853 03854 DO 03855 DO j=1,globalParHeader(-1) 03856 globalParLabelIndex(2,j)=0 ! reset count 03857 END DO 03858 03859 ! read all data files and add all labels to global labels table ---- 03860 03861 IF(mprint /= 0) THEN 03862 WRITE(*,*) 'Read all binary data files:' 03863 END IF 03864 CALL hmpldf(1,'Number of words/record in binary file') 03865 CALL hmpdef(8,0.0,60.0,'not_stored data per record') 03866 ! define read buffer 03867 nc21=ncache/(21*mthrdr) ! split read cache 1 : 10 : 10 for pointers, ints, floats 03868 nwrd=nc21+1 03869 length=nwrd*mthrdr 03870 CALL mpalloc(readBufferPointer,length,'read buffer, pointer') 03871 nwrd=nc21*10+2+ndimbuf 03872 length=nwrd*mthrdr 03873 CALL mpalloc(readBufferDataI,length,'read buffer, integer') 03874 CALL mpalloc(readBufferDataF,length,'read buffer, float') 03875 03876 DO 03877 CALL peread(nr) ! read records 03878 IF (skippedRecords == 0) CALL peprep(0) ! prepare records 03879 IF(nr <= 0) EXIT ! end of data? 03880 END DO 03881 ! release read buffer 03882 CALL mpdealloc(readBufferDataF) 03883 CALL mpdealloc(readBufferDataI) 03884 CALL mpdealloc(readBufferPointer) 03885 IF (skippedRecords == 0) THEN 03886 EXIT 03887 ELSE 03888 WRITE(lunlog,*) 'LOOP1: reading data files again' 03889 END IF 03890 END DO 03891 03892 IF(nhistp /= 0) THEN 03893 CALL hmprnt(1) 03894 CALL hmprnt(8) 03895 END IF 03896 CALL hmpwrt(1) 03897 CALL hmpwrt(8) 03898 CALL upone ! finalize the global label table 03899 ntgb = globalParHeader(-1) ! total number of labels/parameters 03900 WRITE(lunlog,*) 'LOOP1:',ntgb, & 03901 ' is total number NTGB of labels/parameters' 03902 ! histogram number of entries per label ---------------------------- 03903 CALL hmpldf(2,'Number of entries per label') 03904 DO j=1,ntgb 03905 CALL hmplnt(2,globalParLabelIndex(2,j)) 03906 END DO 03907 IF(nhistp /= 0) CALL hmprnt(2) ! print histogram 03908 CALL hmpwrt(2) ! write to his file 03909 03910 ! three subarrays for all global parameters ------------------------ 03911 length=ntgb 03912 CALL mpalloc(globalParameter,length,'global parameters') 03913 globalParameter=0.0D0 03914 CALL mpalloc(globalParPreSigma,length,'pre-sigmas') ! presigmas 03915 globalParPreSigma=0. 03916 CALL mpalloc(globalParStart,length,'global parameters at start') 03917 globalParStart=0. 03918 CALL mpalloc(globalParCopy,length,'copy of global parameters') 03919 03920 DO i=1,lenParameters ! parameter start values 03921 param=listParameters(i)%value 03922 in=inone(listParameters(i)%label) 03923 IF(in /= 0) THEN 03924 globalParameter(in)=param 03925 globalParStart(in)=param 03926 ENDIF 03927 END DO 03928 03929 npresg=0 03930 DO i=1,lenPreSigmas ! pre-sigma values 03931 presg=listPreSigmas(i)%value 03932 in=inone(listPreSigmas(i)%label) 03933 IF(in /= 0) THEN 03934 IF(presg > 0.0) npresg=npresg+1 ! FIXME: check if enough 'entries'? 03935 globalParPreSigma(in)=presg ! insert pre-sigma 0 or > 0 03936 END IF 03937 END DO 03938 WRITE(lunlog,*) 'LOOP1:',npresg,' is number of pre-sigmas' 03939 WRITE(*,*) 'LOOP1:',npresg,' is number of pre-sigmas' 03940 IF(npresg == 0) WRITE(*,*) 'Warning: no pre-sigmas defined' 03941 03942 ! determine flag variable (active) or fixed (inactive) ------------- 03943 03944 indab=0 03945 DO i=1,ntgb 03946 IF(globalParLabelIndex(2,i) >= mreqen.AND.globalParPreSigma(i) >= 0.0) THEN 03947 indab=indab+1 03948 globalParLabelIndex(2,i)=indab ! variable, used in matrix (active) 03949 ELSE 03950 globalParLabelIndex(2,i)=-1 ! fixed, not used in matrix (not active) 03951 END IF 03952 END DO 03953 globalParHeader(-6)=indab ! counted variable 03954 nvgb=indab ! nr of variable parameters 03955 WRITE(lunlog,*) 'LOOP1:',nvgb, ' is number NVGB of variable parameters' 03956 03957 ! translation table of length NVGB of total global indices --------- 03958 length=nvgb 03959 CALL mpalloc(globalParVarToTotal,length,'translation table var -> total') 03960 indab=0 03961 DO i=1,ntgb 03962 IF(globalParLabelIndex(2,i) > 0) THEN 03963 indab=indab+1 03964 globalParVarToTotal(indab)=i 03965 END IF 03966 END DO 03967 03968 ! regularization --------------------------------------------------- 03969 CALL mpalloc(globalParPreWeight,length,'pre-sigmas weights') ! presigma weights 03970 WRITE(*,112) ' Default pre-sigma =',regpre, & 03971 ' (if no individual pre-sigma defined)' 03972 WRITE(*,*) 'Pre-sigma factor is',regula 03973 03974 IF(nregul == 0) THEN 03975 WRITE(*,*) 'No regularization will be done' 03976 ELSE 03977 WRITE(*,*) 'Regularization will be done, using factor',regula 03978 END IF 03979 112 FORMAT(a,e9.2,a) 03980 IF (nvgb <= 0) STOP '... no variable global parameters' 03981 03982 DO ivgbi=1,nvgb ! IVGBI = variable parameter index 03983 itgbi=globalParVarToTotal(ivgbi) ! ITGBI = global parameter index 03984 presg=globalParPreSigma(itgbi) ! get pre-sigma 03985 prewt=0.0 ! pre-weight 03986 IF(presg > 0.0) THEN 03987 prewt=1.0/presg**2 ! 1/presigma^2 03988 ELSE IF(presg == 0.0.AND.regpre > 0.0) THEN 03989 prewt=1.0/regpre**2 ! default 1/presigma^2 03990 END IF 03991 globalParPreWeight(ivgbi)=regula*prewt ! weight = factor / presigma^2 03992 END DO 03993 03994 ! WRITE(*,*) 'GlPa_index GlPa_label array1 array6' 03995 DO i=1,ntgb 03996 itgbl=globalParLabelIndex(1,i) 03997 ivgbi=globalParLabelIndex(2,i) 03998 IF(ivgbi > 0) THEN 03999 ! WRITE(*,111) I,ITGBL,QM(IND1+I),QM(IND6+IVGBI) 04000 ELSE 04001 ! WRITE(*,111) I,ITGBL,QM(IND1+I) 04002 END IF 04003 END DO 04004 ! 111 FORMAT(I5,I10,F10.5,E12.4) 04005 WRITE(*,101) 'NTGB',ntgb,'total number of parameters' 04006 WRITE(*,101) 'NVGB',nvgb,'number of variable parameters' 04007 04008 ! print overview over important numbers ---------------------------- 04009 04010 nrecal=nrec 04011 IF(mprint /= 0) THEN 04012 WRITE(*,*) ' ' 04013 WRITE(*,101) ' NREC',nrec,'number of records' 04014 WRITE(*,101) 'MREQEN',mreqen,'required number of entries' 04015 IF (mreqpe > 1) WRITE(*,101) & 04016 'MREQPE',mreqpe,'required number of pair entries' 04017 IF (msngpe >= 1) WRITE(*,101) & 04018 'MSNGPE',msngpe,'max pair entries single prec. storage' 04019 WRITE(*,101) 'NTGB',ntgb,'total number of parameters' 04020 WRITE(*,101) 'NVGB',nvgb,'number of variable parameters' 04021 IF(mprint > 1) THEN 04022 WRITE(*,*) ' ' 04023 WRITE(*,*) 'Global parameter labels:' 04024 mqi=ntgb 04025 IF(mqi <= 100) THEN 04026 WRITE(*,*) (globalParLabelIndex(2,i),i=1,mqi) 04027 ELSE 04028 WRITE(*,*) (globalParLabelIndex(2,i),i=1,30) 04029 WRITE(*,*) ' ...' 04030 mqi=((mqi-20)/20)*20+1 04031 WRITE(*,*) (globalParLabelIndex(2,i),i=mqi,ntgb) 04032 END IF 04033 END IF 04034 WRITE(*,*) ' ' 04035 WRITE(*,*) ' ' 04036 END IF 04037 WRITE(8,*) ' ' 04038 WRITE(8,101) ' NREC',nrec,'number of records' 04039 WRITE(8,101) 'MREQEN',mreqen,'required number of entries' 04040 04041 WRITE(lunlog,*) 'LOOP1: ending' 04042 WRITE(lunlog,*) ' ' 04043 CALL mend 04044 04045 101 FORMAT(1X,a6,' =',i10,' = ',a) 04046 END SUBROUTINE loop1 04047 04058 04059 SUBROUTINE loop2 04060 USE mpmod 04061 USE mpdalc 04062 04063 IMPLICIT NONE 04064 REAL :: chin2 04065 REAL :: chin3 04066 REAL :: cpr 04067 REAL :: fsum 04068 REAL :: gbc 04069 REAL :: gbu 04070 REAL :: glder 04071 INTEGER :: i 04072 INTEGER :: ia 04073 INTEGER :: ib 04074 INTEGER :: ibuf 04075 INTEGER :: icgb 04076 INTEGER :: iext 04077 INTEGER :: ihis 04078 INTEGER :: ij 04079 INTEGER :: ijn 04080 INTEGER :: inder 04081 INTEGER :: ioff 04082 INTEGER :: iproc 04083 INTEGER :: irecmm 04084 INTEGER :: isfrst 04085 INTEGER :: islast 04086 INTEGER :: ist 04087 INTEGER :: itgbi 04088 INTEGER :: itgbij 04089 INTEGER :: itgbik 04090 INTEGER :: ivgbij 04091 INTEGER :: ivgbik 04092 INTEGER :: j 04093 INTEGER :: ja 04094 INTEGER :: jb 04095 INTEGER :: jcmprs 04096 INTEGER :: jext 04097 INTEGER :: jsp 04098 INTEGER :: k 04099 INTEGER :: kfile 04100 INTEGER :: l 04101 INTEGER :: label 04102 INTEGER :: last 04103 INTEGER :: lu 04104 INTEGER :: lun 04105 INTEGER :: maeqnf 04106 INTEGER :: naeqna 04107 INTEGER :: naeqnf 04108 INTEGER :: naeqng 04109 INTEGER :: nc21 04110 INTEGER :: ncachd 04111 INTEGER :: ncachi 04112 INTEGER :: ncachr 04113 INTEGER :: nda 04114 INTEGER :: ndf 04115 INTEGER :: ndfmax 04116 INTEGER :: nfixed 04117 INTEGER :: nggd 04118 INTEGER :: nggi 04119 INTEGER :: nmatmo 04120 INTEGER :: noff 04121 INTEGER :: nr 04122 INTEGER :: nrecf 04123 INTEGER :: nrecmm 04124 INTEGER :: nst 04125 INTEGER :: nwrd 04126 INTEGER :: inone 04127 REAL :: wgh 04128 REAL :: wolfc3 04129 REAL :: wrec 04130 REAL :: chindl 04131 04132 DOUBLE PRECISION::dstat(3) 04133 INTEGER*8:: noff8 04134 INTEGER(KIND=large):: ndimbi 04135 INTEGER(KIND=large):: ndimsa(4) 04136 INTEGER(KIND=large):: ndgn 04137 INTEGER(KIND=large):: matsiz(2) 04138 INTEGER(KIND=large):: matwords 04139 INTEGER(KIND=large):: length 04140 INTEGER(KIND=large):: rows 04141 INTEGER(KIND=large):: cols 04142 INTEGER(KIND=large), PARAMETER :: two=2 04143 INTEGER:: maxGlobalPar = 0 04144 INTEGER:: maxLocalPar = 0 04145 INTEGER:: maxEquations = 0 04146 04147 INTERFACE ! needed for assumed-shape dummy arguments 04148 SUBROUTINE ndbits(ndims,ncmprs,nsparr,mnpair,ihst,jcmprs) 04149 USE mpdef 04150 INTEGER(KIND=large), DIMENSION(4), INTENT(OUT) :: ndims 04151 INTEGER, DIMENSION(:), INTENT(OUT) :: ncmprs 04152 INTEGER(KIND=large), DIMENSION(:,:), INTENT(OUT) :: nsparr 04153 INTEGER, INTENT(IN) :: mnpair 04154 INTEGER, INTENT(IN) :: ihst 04155 INTEGER, INTENT(IN) :: jcmprs 04156 END SUBROUTINE ndbits 04157 SUBROUTINE spbits(nsparr,nsparc,ncmprs) 04158 USE mpdef 04159 INTEGER(KIND=large), DIMENSION(:,:), INTENT(IN) :: nsparr 04160 INTEGER, DIMENSION(:), INTENT(OUT) :: nsparc 04161 INTEGER, DIMENSION(:), INTENT(IN) :: ncmprs 04162 END SUBROUTINE spbits 04163 END INTERFACE 04164 04165 SAVE 04166 04167 !$ INTEGER :: OMP_GET_THREAD_NUM 04168 04169 isfrst(ibuf)=readBufferPointer(ibuf)+1 04170 islast(ibuf)=readBufferDataI(readBufferPointer(ibuf)) 04171 inder(i)=readBufferDataI(i) 04172 glder(i)=readBufferDataF(i) 04173 ! ... 04174 WRITE(lunlog,*) ' ' 04175 WRITE(lunlog,*) 'LOOP2: starting' 04176 CALL mstart('LOOP2') 04177 04178 ! two subarrays to get the global parameter indices, used in an event 04179 length=nvgb 04180 CALL mpalloc(globalIndexUsage,length,'global index') 04181 CALL mpalloc(backIndexUsage,length,'back index') 04182 backIndexUsage=0 04183 04184 ! constraints - determine number of constraints NCGB 04185 ncgb=0 04186 i=0 04187 last=-1 04188 ! find next constraint header and count nr of constraints 04189 DO WHILE(i < lenConstraints) 04190 i=i+1 04191 label=listConstraints(i)%label 04192 IF(last == 0.AND.label == (-1)) ncgb=ncgb+1 04193 last=label 04194 END DO 04195 WRITE(*,*) 'LOOP2:',ncgb,' constraints' 04196 04197 nagb=nvgb+ncgb ! total number of fit parameters 04198 noff8=int8(nagb)*int8(nagb-1)/2 04199 04200 ! read all data files and add all variable index pairs ------------- 04201 04202 IF(matsto == 2) THEN 04203 IF (mcmprs /= 0) numbit=MAX(numbit,2) ! identify single entries for compression 04204 CALL clbits(nagb,ndimbi,nencdb,numbit) ! get dimension for bit storage, encoding info 04205 END IF 04206 04207 ! reading events===reading events===reading events===reading events= 04208 nrecf =0 ! records with fixed global parameters 04209 naeqng=0 ! count number of equations (with global der.) 04210 naeqnf=0 ! count number of equations ( " , fixed) 04211 naeqna=0 ! all 04212 WRITE(lunlog,*) 'LOOP2: start event reading' 04213 ! monitoring for sparse matrix? 04214 irecmm=0 04215 IF (matsto == 2.AND.matmon /= 0) THEN 04216 nmatmo=0 04217 IF (matmon > 0) THEN 04218 nrecmm=matmon 04219 ELSE 04220 nrecmm=1 04221 END IF 04222 END IF 04223 DO k=1,3 04224 dstat(k)=0.0D0 04225 END DO 04226 ! define read buffer 04227 nc21=ncache/(21*mthrdr) ! split read cache 1 : 10 : 10 for pointers, ints, floats 04228 nwrd=nc21+1 04229 length=nwrd*mthrdr 04230 CALL mpalloc(readBufferPointer,length,'read buffer, pointer') 04231 nwrd=nc21*10+2+ndimbuf 04232 length=nwrd*mthrdr 04233 CALL mpalloc(readBufferDataI,length,'read buffer, integer') 04234 CALL mpalloc(readBufferDataF,length,'read buffer, float') 04235 04236 DO 04237 CALL peread(nr) ! read records 04238 CALL peprep(1) ! prepare records 04239 ioff=0 04240 DO ibuf=1,numReadBuffer ! buffer for current record 04241 nrec=readBufferDataI(isfrst(ibuf)-2) ! record 04242 ! Printout for DEBUG 04243 IF(nrec <= mdebug) THEN 04244 nda=0 04245 kfile=NINT(readBufferDataF(isfrst(ibuf)-1)) ! file 04246 wrec =readBufferDataF(isfrst(ibuf)-2) ! weight 04247 WRITE(*,*) ' ' 04248 WRITE(*,*) 'Record number ',nrec,' from file ',kfile 04249 IF (wgh /= 1.0) WRITE(*,*) ' weight ',wrec 04250 ist=isfrst(ibuf) 04251 nst=islast(ibuf) 04252 DO ! loop over measurements 04253 CALL isjajb(nst,ist,ja,jb,jsp) 04254 IF(ja == 0) EXIT 04255 nda=nda+1 04256 IF(nda > mdebg2) THEN 04257 IF(nda == mdebg2+1) WRITE(*,*) '... and more data' 04258 CYCLE 04259 END IF 04260 WRITE(*,*) ' ' 04261 WRITE(*,*) nda, ' Measured value =',glder(ja),' +- ',glder(jb) 04262 WRITE(*,*) 'Local derivatives:' 04263 WRITE(*,107) (inder(ja+j),glder(ja+j),j=1,jb-ja-1) 04264 107 FORMAT(6(i3,g12.4)) 04265 IF (jb < ist) THEN 04266 WRITE(*,*) 'Global derivatives:' 04267 WRITE(*,108) (globalParLabelIndex(1,inder(jb+j)),inder(jb+j), & 04268 globalParLabelIndex(2,inder(jb+j)),glder(jb+j),j=1,ist-jb) 04269 108 FORMAT(3I11,g12.4) 04270 END IF 04271 IF(nda == 1) THEN 04272 WRITE(*,*) 'total_par_label __label__ var_par_index derivative' 04273 END IF 04274 END DO 04275 WRITE(*,*) ' ' 04276 END IF 04277 04278 nagbn =0 ! count number of global derivatives 04279 nalcn =0 ! count number of local derivatives 04280 naeqn =0 ! count number of equations 04281 maeqnf=naeqnf 04282 ist=isfrst(ibuf) 04283 nst=islast(ibuf) 04284 nwrd=nst-ist+1 04285 DO ! loop over measurements 04286 CALL isjajb(nst,ist,ja,jb,jsp) 04287 IF(ja == 0.AND.jb == 0) EXIT 04288 naeqn=naeqn+1 04289 naeqna=naeqna+1 04290 IF(ja /= 0) THEN 04291 IF (ist > jb) naeqng=naeqng+1 04292 nfixed=0 04293 DO j=1,ist-jb 04294 ij=inder(jb+j) ! index of global parameter 04295 ij=globalParLabelIndex(2,ij) ! change to variable parameter 04296 IF(ij > 0) THEN 04297 ijn=backIndexUsage(ij) ! get index of index 04298 IF(ijn == 0) THEN ! not yet included 04299 nagbn=nagbn+1 ! count 04300 globalIndexUsage(nagbn)=ij ! store variable index 04301 backIndexUsage(ij)=nagbn ! store back index 04302 END IF 04303 ELSE 04304 nfixed=nfixed+1 04305 END IF 04306 END DO 04307 IF (nfixed > 0) naeqnf=naeqnf+1 04308 END IF 04309 04310 IF(ja /= 0.AND.jb /= 0) THEN 04311 DO j=1,jb-ja-1 ! local parameters 04312 ij=inder(ja+j) 04313 nalcn=MAX(nalcn,ij) 04314 END DO 04315 END IF 04316 END DO 04317 04318 ! end-of-event 04319 IF (naeqnf > maeqnf) nrecf=nrecf+1 04320 irecmm=irecmm+1 04321 ! end-of-event-end-of-event-end-of-event-end-of-event-end-of-event-e 04322 04323 maxGlobalPar=MAX(nagbn,maxGlobalPar) ! maximum number of global parameters 04324 maxLocalPar=MAX(nalcn,maxLocalPar) ! maximum number of local parameters 04325 maxEquations=MAX(naeqn,maxEquations) ! maximum number of equations 04326 04327 ! sample statistics for caching 04328 dstat(1)=dstat(1)+dfloat((nwrd+2)*2) ! record size 04329 dstat(2)=dstat(2)+dfloat(nagbn+2) ! indices, 04330 dstat(3)=dstat(3)+dfloat(nagbn*nagbn+nagbn) ! data for MUPDAT 04331 04332 CALL sort1k(globalIndexUsage,nagbn) ! sort global par. 04333 ! overwrite read buffer with lists of global labels 04334 ioff=ioff+1 04335 readBufferPointer(ibuf)=ioff 04336 readBufferDataI(ioff)=ioff+nagbn 04337 DO i=1,nagbn ! reset global index array 04338 iext=globalIndexUsage(i) 04339 backIndexUsage(iext)=0 04340 readBufferDataI(ioff+i)=iext 04341 END DO 04342 ioff=ioff+nagbn 04343 04344 END DO 04345 ioff=0 04346 04347 IF (matsto == 2) THEN 04348 !$OMP PARALLEL & 04349 !$OMP DEFAULT(PRIVATE) & 04350 !$OMP SHARED(numReadBuffer,readBufferPointer,readBufferDataI,MTHRD) 04351 iproc=0 04352 !$ IPROC=OMP_GET_THREAD_NUM() ! thread number 04353 DO ibuf=1,numReadBuffer 04354 ist=isfrst(ibuf) 04355 nst=islast(ibuf) 04356 DO i=ist,nst ! store all combinations 04357 iext=readBufferDataI(i) ! variable global index 04358 !$ IF (MOD(IEXT,MTHRD).EQ.IPROC) THEN ! distinct rows per thread 04359 DO l=ist,i 04360 jext=readBufferDataI(l) 04361 CALL inbits(iext,jext,1) ! save space 04362 END DO 04363 !$ ENDIF 04364 END DO 04365 END DO 04366 !$OMP END PARALLEL 04367 ! monitoring 04368 IF (matmon /= 0.AND. & 04369 (irecmm >= nrecmm.OR.irecmm == mxrec)) THEN 04370 IF (nmatmo == 0) THEN 04371 WRITE(*,*) 04372 WRITE(*,*) 'Monitoring of sparse matrix construction' 04373 WRITE(*,*) ' records ........ off-diagonal elements ', & 04374 '....... compression memory' 04375 WRITE(*,*) ' non-zero used(double) used', & 04376 '(float) [%] [GB]' 04377 END IF 04378 nmatmo=nmatmo+1 04379 jcmprs=MAX(mcmprs,msngpe) 04380 CALL ckbits(ndimsa,mreqpe,jcmprs) 04381 gbc=4.0E-9*FLOAT(ndimsa(2)+2*ndimsa(3)+ndimsa(4)) ! GB compressed 04382 gbu=1.2E-8*FLOAT(ndimsa(3)+ndimsa(4)) ! GB uncompressed 04383 cpr=100.0*gbc/gbu 04384 WRITE(*,1177) irecmm,ndimsa(1),ndimsa(3),ndimsa(4),cpr,gbc 04385 1177 FORMAT(i9,3I13,f10.2,f11.6) 04386 DO WHILE(irecmm >= nrecmm) 04387 IF (matmon > 0) THEN 04388 nrecmm=nrecmm+matmon 04389 ELSE 04390 nrecmm=nrecmm*2 04391 END IF 04392 END DO 04393 END IF 04394 04395 END IF 04396 04397 IF (nr <= 0) EXIT ! next block of events ? 04398 END DO 04399 ! release read buffer 04400 CALL mpdealloc(readBufferDataF) 04401 CALL mpdealloc(readBufferDataI) 04402 CALL mpdealloc(readBufferPointer) 04403 04404 WRITE(lunlog,*) 'LOOP2: event reading ended - end of data' 04405 DO k=1,3 04406 dstat(k)=dstat(k)/dfloat(nrec) 04407 END DO 04408 ! end=of=data=end=of=data=end=of=data=end=of=data=end=of=data=end=of 04409 04410 IF(matsto == 2) THEN 04411 04412 ! constraints and index pairs with Lagrange multiplier 04413 04414 04415 ! constraints - determine number of constraints NCGB and index-pairs 04416 ! Lagrange multiplier and global parameters 04417 04418 04419 i=0 04420 icgb=0 04421 last=-1 04422 ! find next constraint header 04423 DO WHILE(i < lenConstraints) 04424 i=i+1 04425 label=listConstraints(i)%label 04426 IF(last == 0.AND.label == (-1)) icgb=icgb+1 04427 IF(label > 0) THEN 04428 itgbi=inone(label) 04429 ij=globalParLabelIndex(2,itgbi) ! change to variable parameter 04430 IF(ij > 0) THEN 04431 CALL inbits(nvgb+icgb,ij,mreqpe) 04432 END IF 04433 END IF 04434 last=label 04435 END DO 04436 04437 ! measurements - determine index-pairs 04438 04439 04440 i=1 04441 DO WHILE (i <= lenMeasurements) 04442 i=i+2 04443 ! loop over label/factor pairs 04444 ia=i 04445 DO 04446 i=i+1 04447 IF(i > lenMeasurements) EXIT 04448 IF(listMeasurements(i)%label == 0) EXIT 04449 END DO 04450 ib=i-1 04451 04452 DO j=ia,ib 04453 itgbij=inone(listMeasurements(j)%label) ! total parameter index 04454 ! first index 04455 ivgbij=0 04456 IF(itgbij /= 0) ivgbij=globalParLabelIndex(2,itgbij) ! variable-parameter index 04457 DO k=ia,j 04458 itgbik=inone(listMeasurements(k)%label) ! total parameter index 04459 ! second index 04460 ivgbik=0 04461 IF(itgbik /= 0) ivgbik=globalParLabelIndex(2,itgbik) ! variable-parameter index 04462 IF(ivgbij > 0.AND.ivgbik > 0) THEN 04463 CALL inbits(ivgbij,ivgbik,mreqpe) 04464 IF (mprint > 1) WRITE(*,*) 'add index pair ',ivgbij,ivgbik 04465 END IF 04466 END DO 04467 END DO 04468 04469 END DO 04470 END IF 04471 04472 ! print numbers ---------------------------------------------------- 04473 04474 IF (nagb >= 65536) THEN 04475 noff=INT(noff8/1000) 04476 ELSE 04477 noff=INT(noff8) 04478 END IF 04479 ndgn=0 04480 nspc=1 ! number of precision types (double, single) for matrix storage 04481 matwords=0 04482 IF(matsto == 2) THEN 04483 ihis=0 04484 IF (mhispe > 0) THEN 04485 ihis=15 04486 CALL hmpdef(ihis,0.0,FLOAT(mhispe), 'NDBITS: #off-diagonal elements') 04487 END IF 04488 jcmprs=MAX(mcmprs,msngpe) 04489 IF (jcmprs > 0.AND.numbit > 1) nspc=2 ! mixed precision storage 04490 length=nagb*nspc 04491 CALL mpalloc(sparseMatrixCompression,length,'INBITS: row compression') 04492 sparseMatrixCompression=0 04493 length=(nagb+1)*nspc 04494 CALL mpalloc(sparseMatrixOffsets,two,length, 'sparse matrix row offsets') 04495 CALL ndbits(ndimsa,sparseMatrixCompression,sparseMatrixOffsets, & 04496 mreqpe,ihis,jcmprs) 04497 ndgn=ndimsa(3)+ndimsa(4) ! actual number of off-diagonal elements 04498 matwords=ndimsa(2)+length ! size of sparsity structure 04499 04500 IF (mhispe > 0) THEN 04501 IF (nhistp /= 0) CALL hmprnt(ihis) 04502 CALL hmpwrt(ihis) 04503 END IF 04504 END IF 04505 04506 nagbn=maxGlobalPar ! max number of global parameters in one event 04507 nalcn=maxLocalPar ! max number of local parameters in one event 04508 naeqn=maxEquations ! max number of equations in one event 04509 CALL mpdealloc(globalIndexUsage) 04510 CALL mpdealloc(backIndexUsage) 04511 ! matrices for event matrices 04512 ! split up cache 04513 IF (fcache(2) == 0.0) THEN ! from data (DSTAT) 04514 fcache(1)=SNGL(dstat(1))*fcache(1) ! leave some part free for fluctuations 04515 fcache(2)=SNGL(dstat(2)) 04516 fcache(3)=SNGL(dstat(3)) 04517 END IF 04518 fsum=fcache(1)+fcache(2)+fcache(3) 04519 DO k=1,3 04520 fcache(k)=fcache(k)/fsum 04521 END DO 04522 ncachr=IFIX(FLOAT(ncache)*fcache(1)+0.5) ! read cache 04523 ! define read buffer 04524 nc21=ncachr/(21*mthrdr) ! split read cache 1 : 10 : 10 for pointers, ints, floats 04525 nwrd=nc21+1 04526 length=nwrd*mthrdr 04527 CALL mpalloc(readBufferPointer,length,'read buffer, pointer') 04528 nwrd=nc21*10+2+ndimbuf 04529 length=nwrd*mthrdr 04530 CALL mpalloc(readBufferDataI,length,'read buffer, integer') 04531 CALL mpalloc(readBufferDataF,length,'read buffer, float') 04532 04533 ncachi=IFIX(FLOAT(ncache)*fcache(2)+0.5) ! index cache 04534 ncachd=ncache-ncachr-ncachi ! data cache 04535 nggd=(nagbn*nagbn+nagbn)/2+ncachd/(2*mthrd) ! number of double 04536 nggi=2+nagbn+ncachi/mthrd ! number of ints 04537 length=nagbn*mthrd 04538 CALL mpalloc(globalIndexUsage,length, 'global parameters (dim =max/event)') 04539 length=nvgb*mthrd 04540 CALL mpalloc(backIndexUsage,length,'global variable-index array') 04541 backIndexUsage=0 04542 length=nagbn*nalcn 04543 CALL mpalloc(localGlobalMatrix,length,'local/global matrix') 04544 length=nggd*mthrd 04545 CALL mpalloc(writeBufferUpdates,length,'symmetric update matrices') 04546 writeBufferHeader(-1)=nggd ! number of words per thread 04547 writeBufferHeader(-2)=(nagbn*nagbn+nagbn)/2 ! min free (double) words 04548 length=nggi*mthrd 04549 CALL mpalloc(writeBufferIndices,length,'symmetric update matrix indices') 04550 rows=7; cols=mthrd 04551 CALL mpalloc(writeBufferInfo,rows,cols,'write buffer status (I)') 04552 rows=2; cols=mthrd 04553 CALL mpalloc(writeBufferData,rows,cols,'write buffer status (F)') 04554 writeBufferHeader(1)=nggi ! number of words per thread 04555 writeBufferHeader(2)=nagbn+2 ! min free words 04556 04557 ! print all relevant dimension parameters 04558 04559 DO lu=6,8,2 ! unit 6 and 8 04560 04561 WRITE(*,*) ' ' 04562 WRITE(lu,101) 'NTGB',ntgb,'total number of parameters' 04563 WRITE(lu,102) '(all parameters, appearing in binary files)' 04564 WRITE(lu,101) 'NVGB',nvgb,'number of variable parameters' 04565 WRITE(lu,102) '(appearing in fit matrix/vectors)' 04566 WRITE(lu,101) 'NAGB',nagb,'number of fit parameters' 04567 WRITE(lu,102) '(including Lagrange multiplier or reduced)' 04568 WRITE(lu,101) 'MBANDW',mbandw,'band width of band matrix' 04569 WRITE(lu,102) '(if =0, no band matrix)' 04570 IF (nagb >= 65536) THEN 04571 WRITE(lu,101) 'NOFF/K',noff,'max number of off-diagonal elements' 04572 ELSE 04573 WRITE(lu,101) 'NOFF',noff,'max number of off-diagonal elements' 04574 END IF 04575 IF(ndgn /= 0) & 04576 WRITE(lu,101) 'NDGN',ndgn,'actual number of off-diagonal elements' 04577 WRITE(lu,101) 'NCGB',ncgb,'number of constraints' 04578 WRITE(lu,101) 'NAGBN',nagbn,'max number of global parameters in an event' 04579 WRITE(lu,101) 'NALCN',nalcn,'max number of local parameters in an event' 04580 WRITE(lu,101) 'NAEQN',naeqn,'max number of equations in an event' 04581 IF (mprint > 1) THEN 04582 WRITE(lu,101) 'NAEQNA',naeqna,'number of equations' 04583 WRITE(lu,101) 'NAEQNG',naeqng, & 04584 'number of equations with global derivatives' 04585 WRITE(lu,101) 'NAEQNF',naeqnf, & 04586 'number of equations with fixed global derivatives' 04587 WRITE(lu,101) 'NRECF',nrecf, & 04588 'number of records with fixed global derivatives' 04589 END IF 04590 IF (ncache > 0) THEN 04591 WRITE(lu,101) 'NCACHE',ncache,'number of words for caching' 04592 WRITE(lu,111) (fcache(k)*100.0,k=1,3) 04593 111 FORMAT(22X,'cache splitting ',3(f6.1,' %')) 04594 END IF 04595 WRITE(lu,*) ' ' 04596 IF(matsto == 2) THEN 04597 ! FRUSE=FLOAT(MQ(INDF-1))/FLOAT(NOFF) 04598 ! WRITE(LU,123) 100.0*FRUSE 04599 END IF 04600 ! 123 FORMAT(' Fraction of non-zero off-diagonal elements is', 04601 ! + F6.1,' %') 04602 04603 ! print information about methods and matrix storage modes 04604 04605 WRITE(lu,*) ' ' 04606 WRITE(lu,*) 'Solution method and matrix-storage mode:' 04607 IF(metsol == 1) THEN 04608 WRITE(lu,*) ' METSOL = 1: matrix inversion' 04609 ELSE IF(metsol == 2) THEN 04610 WRITE(lu,*) ' METSOL = 2: diagonalization' 04611 ELSE IF(metsol == 3) THEN 04612 !GF WRITE(LU,*) ' METSOL = 3: MINRES' 04613 WRITE(lu,*) ' METSOL = 3: MINRES (rtol', mrestl,')' 04614 ELSE IF(metsol == 4) THEN 04615 WRITE(lu,*) ' METSOL = 4: GMRES' 04616 END IF 04617 WRITE(lu,*) ' with',mitera,' iterations' 04618 IF(matsto == 1) THEN 04619 WRITE(lu,*) ' MATSTO = 1: symmetric matrix, ', '(n*n+n)/2 elements' 04620 ELSE IF(matsto == 2) THEN 04621 WRITE(lu,*) ' MATSTO = 2: sparse matrix' 04622 END IF 04623 ! END IF 04624 IF(dflim /= 0.0) THEN 04625 WRITE(lu,103) 'Convergence assumed, if expected dF <',dflim 04626 END IF 04627 04628 END DO ! print loop 04629 04630 ! Wolfe conditions 04631 04632 IF(0.0 < wolfc1.AND.wolfc1 < wolfc2.AND.wolfc2 < 1.0) GO TO 32 04633 IF(wolfc1 == 0.0) wolfc1=1.0E-4 04634 IF(wolfc2 == 0.0) wolfc2=0.9 04635 IF(0.0 < wolfc1.AND.wolfc1 < wolfc2.AND.wolfc2 < 1.0) GO TO 32 04636 IF(wolfc1 <= 0.0) wolfc1=1.0E-4 04637 IF(wolfc2 >= 1.0) wolfc2=0.9 04638 IF(wolfc1 > wolfc2) THEN ! exchange 04639 wolfc3=wolfc1 04640 wolfc1=wolfc2 04641 wolfc2=wolfc3 04642 ELSE 04643 wolfc1=1.0E-4 04644 wolfc2=0.9 04645 END IF 04646 WRITE(*,105) wolfc1,wolfc2 04647 WRITE(lun,105) wolfc1,wolfc2 04648 105 FORMAT(' Constants C1, C2 for Wolfe conditions:',g12.4,', ',g12.4) 04649 04650 ! prepare matrix and gradient storage ------------------------------ 04651 !32 CONTINUE 04652 32 matsiz(1)=int8(nagb)*int8(nagb+1)/2 ! number of words for double precision storage 'j' 04653 matsiz(2)=0 ! number of words for single precision storage '.' 04654 IF(matsto == 2) THEN ! sparse matrix 04655 matsiz(1)=ndimsa(3)+nagb 04656 matsiz(2)=ndimsa(4) 04657 CALL mpalloc(sparseMatrixColumns,ndimsa(2),'sparse matrix column list') 04658 CALL spbits(sparseMatrixOffsets,sparseMatrixColumns,sparseMatrixCompression) 04659 END IF 04660 matwords=matwords+matsiz(1)*2+matsiz(2) ! #words for matrix storage 04661 04662 CALL feasma ! prepare constraint matrices 04663 04664 CALL vmprep(matsiz) ! prepare matrix and gradient storage 04665 WRITE(*,*) ' ' 04666 IF (matwords < 250000) THEN 04667 WRITE(*,*) 'Size of global matrix: < 1 MB' 04668 ELSE 04669 WRITE(*,*) 'Size of global matrix:',INT(FLOAT(matwords)*4.0E-6),' MB' 04670 ENDIF 04671 ! print chi^2 cut tables 04672 04673 ndfmax=naeqn-1 04674 WRITE(lunlog,*) ' ' 04675 WRITE(lunlog,*) ' Cut values of Chi^2/Ndf and Chi2,' 04676 WRITE(lunlog,*) ' corresponding to 2 and 3 standard deviations' 04677 WRITE(lunlog,*) ' Ndf Chi^2/Ndf(2) Chi^2(2) ', & 04678 ' Chi^2/Ndf(3) Chi^2(3)' 04679 ndf=0 04680 DO 04681 IF(ndf > naeqn) EXIT 04682 IF(ndf < 10) THEN 04683 ndf=ndf+1 04684 ELSE IF(ndf < 20) THEN 04685 ndf=ndf+2 04686 ELSE IF(ndf < 100) THEN 04687 ndf=ndf+5 04688 ELSE IF(ndf < 200) THEN 04689 ndf=ndf+10 04690 ELSE 04691 EXIT 04692 END IF 04693 chin2=chindl(2,ndf) 04694 chin3=chindl(3,ndf) 04695 WRITE(lunlog,106) ndf,chin2,chin2*FLOAT(ndf),chin3, chin3*FLOAT(ndf) 04696 END DO 04697 04698 WRITE(lunlog,*) 'LOOP2: ending' 04699 WRITE(lunlog,*) ' ' 04700 CALL mend 04701 101 FORMAT(1X,a6,' =',i10,' = ',a) 04702 102 FORMAT(22X,a) 04703 103 FORMAT(1X,a,g12.4) 04704 106 FORMAT(i6,2(3X,f9.3,f12.1,3X)) 04705 END SUBROUTINE loop2 04706 04710 04711 SUBROUTINE vmprep(msize) 04712 USE mpmod 04713 USE mpdalc 04714 04715 IMPLICIT NONE 04716 INTEGER :: i 04717 INTEGER :: nmats 04718 ! 04719 INTEGER(KIND=large), INTENT(IN) :: msize(2) 04720 04721 INTEGER(KIND=large) :: length 04722 SAVE 04723 ! ... 04724 ! Vector/matrix storage 04725 length=nagb*mthrd 04726 CALL mpalloc(globalVector,length,'rhs vector') ! double precision vector 04727 lenGlobalVec=nagb 04728 length=naeqn 04729 CALL mpalloc(localCorrections,length,'residual vector of one record') 04730 length=nalcn*nalcn 04731 CALL mpalloc(aux,length,' local fit scratch array: aux') 04732 CALL mpalloc(vbnd,length,' local fit scratch array: vbnd') 04733 CALL mpalloc(vbdr,length,' local fit scratch array: vbdr') 04734 length=((nalcn+1)*nalcn)/2 04735 CALL mpalloc(clmat,length,' local fit matrix: clmat') 04736 CALL mpalloc(vbk,length,' local fit scratch array: vbk') 04737 length=nalcn 04738 CALL mpalloc(blvec,length,' local fit vector: blvec') 04739 CALL mpalloc(vzru,length,' local fit scratch array: vzru') 04740 CALL mpalloc(scdiag,length,' local fit scratch array: scdiag') 04741 CALL mpalloc(scflag,length,' local fit scratch array: scflag') 04742 CALL mpalloc(ibandh,length,' local fit band width hist.: ibandh') 04743 04744 CALL mpalloc(globalMatD,msize(1),'global matrix (D)' ) 04745 CALL mpalloc(globalMatF,msize(2),'global matrix (F)') 04746 04747 IF(metsol >= 3) THEN ! GMRES/MINRES algorithms 04748 ! array space is: 04749 ! variable-width band matrix or diagonal matrix for parameters 04750 ! followed by rectangular matrix for constraints 04751 ! followed by symmetric matrix for constraints 04752 IF(mbandw > 0) THEN ! variable-width band matrix 04753 length=nagb 04754 CALL mpalloc(indPreCond,length,'pointer-array variable-band matrix') 04755 DO i=1,MIN(mbandw,nvgb) 04756 indPreCond(i)=(i*i+i)/2 ! increasing number 04757 END DO 04758 DO i=MIN(mbandw,nvgb)+1,nvgb 04759 indPreCond(i)=indPreCond(i-1)+mbandw ! fixed band width 04760 END DO 04761 DO i=nvgb+1,nagb ! reset 04762 indPreCond(i)=0 04763 END DO 04764 length=indPreCond(nvgb)+ncgb*nvgb+(ncgb*ncgb+ncgb)/2 04765 CALL mpalloc(matPreCond,length,'variable-band matrix') 04766 ELSE ! default preconditioner 04767 length=nvgb+ncgb*nvgb+(ncgb*ncgb+ncgb)/2 04768 CALL mpalloc(matPreCond,length,'default preconditioner matrix') 04769 END IF 04770 END IF 04771 04772 04773 length=nagb 04774 CALL mpalloc(globalCorrections,length,'corrections') ! double prec corrections 04775 04776 CALL mpalloc(workspaceLinesearch,length,'auxiliary array (D2)') ! double aux 2 04777 CALL mpalloc(workspaceI, length,'auxiliary array (I)') ! int aux 1 04778 CALL mpalloc(workspaceMinres,length,'auxiliary array (D7)') ! double aux 7 04779 04780 IF(metsol == 1) THEN 04781 CALL mpalloc(workspaceD,length,'auxiliary array (D1)') ! double aux 1 04782 ! CALL MEGARR('t D',2*NAGB,'auxiliary array') ! double aux 8 04783 END IF 04784 04785 IF(metsol == 2) THEN 04786 CALL mpalloc(workspaceD,length,'auxiliary array (D1)') ! double aux 1 04787 CALL mpalloc(workspaceDiagonalization,length,'auxiliary array (D3)') ! double aux 3 04788 CALL mpalloc(workspaceEigenValues,length,'auxiliary array (D6)') ! double aux 6 04789 length=nagb*nagb 04790 CALL mpalloc(workspaceEigenVectors,length,'(rotation) matrix U') ! rotation matrix 04791 END IF 04792 04793 IF(metsol >= 3) THEN 04794 nmats=6 04795 length=nmats*nagb 04796 CALL mpalloc(workspaceD,length,'auxiliary array (D1)') ! double aux 1 04797 END IF 04798 04799 END SUBROUTINE vmprep 04800 04804 04805 SUBROUTINE minver 04806 USE mpmod 04807 04808 IMPLICIT NONE 04809 INTEGER :: lun 04810 INTEGER :: ndefec 04811 INTEGER :: nrank 04812 04813 SAVE 04814 ! ... 04815 lun=lunlog ! log file 04816 IF(lunlog == 0) lunlog=6 04817 04818 ! WRITE(*,*) 'MINVER ICALCM=',ICALCM 04819 IF(icalcm == 1) THEN 04820 CALL sqminl(globalMatD, globalCorrections,nagb,nrank, & 04821 workspaceD,workspaceI) 04822 ndefec=nagb-nrank ! rank defect 04823 IF(ndefec /= 0) THEN 04824 WRITE(*,*) 'The rank defect of the symmetric',nagb, & 04825 '-by-',nagb,' matrix is ',ndefec,' (should be zero).' 04826 WRITE(lun,*) 'The rank defect of the symmetric',nagb, & 04827 '-by-',nagb,' matrix is ',ndefec,' (should be zero).' 04828 IF (iforce == 0) THEN 04829 isubit=1 04830 WRITE(*,*) ' --> enforcing SUBITO mode' 04831 WRITE(lun,*) ' --> enforcing SUBITO mode' 04832 END IF 04833 ELSE 04834 WRITE(lun,*) 'No rank defect of the symmetric matrix' 04835 END IF 04836 04837 ELSE ! multiply gradient by inverse matrix 04838 CALL dbsvxl(globalMatD,globalVector,globalCorrections,nagb) 04839 END IF 04840 END SUBROUTINE minver 04841 04843 SUBROUTINE mdiags 04844 USE mpmod 04845 04846 IMPLICIT NONE 04847 REAL :: evalue 04848 INTEGER :: i 04849 INTEGER :: iast 04850 INTEGER :: idia 04851 INTEGER :: imin 04852 INTEGER :: lun 04853 INTEGER :: nmax 04854 INTEGER :: nmin 04855 INTEGER :: ntop 04856 INTEGER :: nvar 04857 ! 04858 INTEGER(KIND=large) :: ii 04859 SAVE 04860 ! ... 04861 04862 lun=lunlog ! log file 04863 IF(lunlog == 0) lun=6 04864 04865 IF(icalcm == 1) THEN 04866 nvar=nagb 04867 DO i=1,nvar ! used in FEASIB 04868 ii=i 04869 workspaceD(i)=globalMatD((ii*ii+ii)/2) ! save diagonal elements 04870 END DO 04871 04872 ! eigenvalues eigenvectors symm_input 04873 CALL devrot(nvar,workspaceEigenValues,workspaceEigenVectors,globalMatD, & 04874 workspaceDiagonalization,workspaceI) 04875 04876 ! histogram of positive eigenvalues 04877 04878 nmax=INT(1.0+LOG10(SNGL(workspaceEigenValues(1)))) ! > log of largest eigenvalue 04879 imin=1 04880 DO i=nagb,1,-1 04881 IF(workspaceEigenValues(i) > 0.0D0) THEN 04882 imin=i ! index of smallest pos. eigenvalue 04883 EXIT 04884 END IF 04885 END DO 04886 nmin=INT(LOG10(SNGL(workspaceEigenValues(imin)))) ! log of smallest pos. eigenvalue 04887 ntop=nmin+6 04888 DO WHILE(ntop < nmax) 04889 ntop=ntop+3 04890 END DO 04891 04892 CALL hmpdef(7,FLOAT(nmin),FLOAT(ntop), 'log10 of positive eigenvalues') 04893 DO idia=1,nagb 04894 IF(workspaceEigenValues(idia) > 0.0D0) THEN ! positive 04895 evalue=LOG10(SNGL(workspaceEigenValues(idia))) 04896 CALL hmpent(7,evalue) 04897 END IF 04898 END DO 04899 IF(nhistp /= 0) CALL hmprnt(7) 04900 CALL hmpwrt(7) 04901 04902 iast=MAX(1,imin-60) 04903 CALL gmpdef(3,2,'low-value end of eigenvalues') 04904 DO i=iast,nagb 04905 evalue=SNGL(workspaceEigenValues(i)) 04906 CALL gmpxy(3,FLOAT(i),evalue) 04907 END DO 04908 IF(nhistp /= 0) CALL gmprnt(3) 04909 CALL gmpwrt(3) 04910 04911 DO i=1,nvar 04912 workspaceDiagonalization(i)=0.0D0 04913 IF(workspaceEigenValues(i) /= 0.0D0) THEN 04914 workspaceDiagonalization(i)=MAX(0.0D0,LOG10(ABS(workspaceEigenValues(i)))+3.0D0) 04915 IF(workspaceEigenValues(i) < 0.0D0) workspaceDiagonalization(i)=-workspaceDiagonalization(i) 04916 END IF 04917 END DO 04918 WRITE(lun,*) ' ' 04919 WRITE(lun,*) 'The first (largest) eigenvalues ...' 04920 WRITE(lun,102) (workspaceEigenValues(i),i=1,MIN(20,nagb)) 04921 WRITE(lun,*) ' ' 04922 WRITE(lun,*) 'The last eigenvalues ... up to',nvgb 04923 WRITE(lun,102) (workspaceEigenValues(i),i=MAX(1,nvgb-19),nvgb) 04924 WRITE(lun,*) ' ' 04925 WRITE(lun,*) 'The eigenvalues from',nvgb+1,' to',nagb 04926 WRITE(lun,102) (workspaceEigenValues(i),i=nvgb+1,nagb) 04927 WRITE(lun,*) ' ' 04928 WRITE(lun,*) 'Log10 + 3 of ',nagb,' eigenvalues in decreasing', ' order' 04929 WRITE(lun,*) '(for Eigenvalue < 0.001 the value 0.0 is shown)' 04930 WRITE(lun,101) (workspaceDiagonalization(i),i=1,nagb) 04931 IF(workspaceDiagonalization(nvar) < 0) WRITE(lun,*) 'Negative values are ', & 04932 'printed for negative eigenvalues' 04933 CALL devsig(nagb,workspaceEigenValues,workspaceEigenVectors,globalVector,workspaceDiagonalization) 04934 WRITE(lun,*) ' ' 04935 WRITE(lun,*) nvgb,' significances: insignificant if ', & 04936 'compatible with N(0,1)' 04937 WRITE(lun,101) (workspaceDiagonalization(i),i=1,nvgb) 04938 04939 04940 101 FORMAT(10F7.1) 04941 102 FORMAT(5E14.6) 04942 04943 END IF 04944 04945 ! solution --------------------------------------------------------- 04946 04947 ! eigenvalues eigenvectors 04948 CALL devsol(nvar,workspaceEigenValues,workspaceEigenVectors,globalVector,globalCorrections,workspaceDiagonalization) 04949 RETURN 04950 END SUBROUTINE mdiags 04951 04953 SUBROUTINE zdiags 04954 USE mpmod 04955 04956 IMPLICIT NONE 04957 ! eigenvalue eigenvectors cov.matrix 04958 CALL devinv(nagb,workspaceEigenValues,workspaceEigenVectors,globalMatD) ! inv 04959 04960 END SUBROUTINE zdiags 04961 04967 04968 SUBROUTINE mminrs 04969 USE mpmod 04970 04971 IMPLICIT NONE 04972 INTEGER :: istop 04973 INTEGER :: itn 04974 INTEGER :: itnlim 04975 INTEGER :: lun 04976 INTEGER :: nout 04977 INTEGER :: nrkd 04978 INTEGER :: nrkd2 04979 04980 DOUBLE PRECISION :: shift 04981 DOUBLE PRECISION :: rtol 04982 DOUBLE PRECISION :: anorm 04983 DOUBLE PRECISION :: acond 04984 DOUBLE PRECISION :: rnorm 04985 DOUBLE PRECISION :: ynorm 04986 LOGICAL :: checka 04987 EXTERNAL avprod, mvsolv, mcsolv 04988 SAVE 04989 ! ... 04990 lun=lunlog ! log file 04991 IF(lunlog == 0) lun=6 04992 04993 nout=lun 04994 itnlim=2000 ! iteration limit 04995 shift =0.0D0 ! not used 04996 rtol = mrestl ! from steering 04997 checka=.FALSE. 04998 04999 ! WRITE(*,*) 'METSOL MATSTO',METSOL,MATSTO 05000 05001 ! WRITE(*,*) 'Input vector' 05002 ! WRITE(*,*) (DQ(IGVEC/2+I),I=1,NAGB) 05003 05004 IF(mbandw == 0) THEN ! default preconditioner 05005 IF(icalcm == 1) THEN 05006 ! WRITE(LUN,*) 'MMINRS: PRECON started' 05007 CALL precon(ncgb,nvgb,matPreCond,matPreCond, matPreCond(1+nvgb), & 05008 matPreCond(1+nvgb+ncgb*nvgb)) 05009 ! WRITE(LUN,*) 'MMINRS: PRECON ended' 05010 END IF 05011 05012 ! WRITE(*,*) 'NAGB MBANDW ',NAGB,MBANDW 05013 ! WRITE(*,*) IGVEC,IDUX1,ISOLV,IDUX7 05014 05015 CALL minres(nagb,globalVector, & 05016 workspaceD(1), workspaceD(1+nagb), workspaceD(1+2*nagb), & 05017 workspaceD(1+3*nagb),workspaceD(1+4*nagb),workspaceD(1+5*nagb), & 05018 globalCorrections,workspaceMinres, avprod, mcsolv, checka ,.TRUE. , shift, & 05019 nout , itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm) 05020 ! WRITE(*,*) 'MINRES ended' 05021 ! + DQ(IDUX1/2+1),DQ(IDUX8/2+1),DQ(IDUX3/2+1), 05022 ! + DQ(IDUX4/2+1),DQ(IDUX5/2+1),DQ(IDUX6/2+1), 05023 ELSE IF(mbandw > 0) THEN ! band matrix preconditioner 05024 IF(icalcm == 1) THEN 05025 WRITE(lun,*) 'MMINRS: EQUDEC started' 05026 CALL equdec(nvgb,ncgb,matPreCond,indPreCond,nrkd,nrkd2) 05027 WRITE(lun,*) 'MMINRS: EQUDEC ended' 05028 END IF 05029 CALL minres(nagb,globalVector, & 05030 workspaceD(1), workspaceD(1+nagb), workspaceD(1+2*nagb), & 05031 workspaceD(1+3*nagb),workspaceD(1+4*nagb),workspaceD(1+5*nagb), & 05032 globalCorrections,workspaceMinres, avprod, mvsolv, checka ,.TRUE. , shift, & 05033 nout, itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm) 05034 ELSE 05035 ! IF(ICALCM.NE.0) WRITE(LUN,*) 'MMINRS: no preconditioning !!!' 05036 CALL minres(nagb,globalVector, & 05037 workspaceD(1), workspaceD(1+nagb), workspaceD(1+2*nagb), & 05038 workspaceD(1+3*nagb),workspaceD(1+4*nagb),workspaceD(1+5*nagb), & 05039 globalCorrections,workspaceMinres, avprod, mvsolv, checka ,.FALSE., shift, & 05040 nout, itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm) 05041 END IF 05042 iitera=itn 05043 istopa=istop 05044 mnrsit=mnrsit+itn 05045 05046 IF (istopa == 0) PRINT *, 'MINRES: istop=0, exact solution x=0.' 05047 ! WRITE(*,*) 'Solution vector' 05048 ! WRITE(*,*) (DQ(ISOLV/2+I),I=1,NAGB) 05049 05050 END SUBROUTINE mminrs 05051 05059 05060 SUBROUTINE mcsolv(n,x,y) ! solve M*y = x 05061 USE mpmod 05062 05063 IMPLICIT NONE 05064 INTEGER,INTENT(IN) :: n 05065 DOUBLE PRECISION, INTENT(IN) :: x(n) 05066 DOUBLE PRECISION, INTENT(OUT) :: y(n) 05067 SAVE 05068 ! ... 05069 CALL presol(ncgb,nvgb,matPreCond,matPreCond(1+nvgb),matPreCond(1+nvgb+ncgb*nvgb),y,x) 05070 END SUBROUTINE mcsolv 05071 05079 05080 SUBROUTINE mvsolv(n,x,y) ! solve M*y = x 05081 USE mpmod 05082 05083 IMPLICIT NONE 05084 REAL :: factr 05085 INTEGER :: i 05086 INTEGER :: icgb 05087 INTEGER :: itgbi 05088 INTEGER :: ivgb 05089 INTEGER :: label 05090 INTEGER :: last 05091 INTEGER :: inone 05092 05093 INTEGER, INTENT(IN) :: n 05094 DOUBLE PRECISION, INTENT(IN) :: x(n) 05095 DOUBLE PRECISION, INTENT(OUT) :: y(n) 05096 SAVE 05097 ! ... 05098 DO i=1,n 05099 y(i)=x(i) ! copy to output vector 05100 END DO 05101 05102 DO icgb=1,ncgb ! copy previous lambda values 05103 ! Y(NVGB+ICGB)=DQ(ISOLV/2+NVGB+ICGB) 05104 END DO 05105 05106 i=0 05107 icgb=0 05108 last=-1 05109 DO WHILE(i < lenConstraints) 05110 i=i+1 05111 label=listConstraints(i)%label 05112 IF(last == 0.AND.label == (-1)) icgb=icgb+1 ! next constraint header 05113 IF(label > 0) THEN 05114 factr=listConstraints(i)%value 05115 itgbi=inone(label) ! -> ITGBI= index of parameter label 05116 ivgb =globalParLabelIndex(2,itgbi) ! -> variable-parameter index 05117 ! - A^T * lambda 05118 IF(ivgb /= 0) THEN 05119 ! Y(IVGB)=Y(IVGB) - DQ(ISOLV/2+NVGB+ICGB)*FACTR 05120 END IF 05121 END IF 05122 last=label 05123 END DO 05124 05125 ! CALL VABSLV(NVGB,DQ(INDV/2+1),MQ(INDU+1),Y) ! solve for Y 05126 ! CALL VABSLV(NAGB,DQ(INDV/2+1),MQ(INDU+1),Y) ! solve for Y 05127 05128 CALL equslv(nvgb,ncgb,matPreCond,indPreCond,y) 05129 END SUBROUTINE mvsolv 05130 05131 05132 05133 !*********************************************************************** 05134 05147 05148 SUBROUTINE xloopn ! 05149 USE mpmod 05150 05151 IMPLICIT NONE 05152 REAL :: catio 05153 REAL :: concu2 05154 REAL :: concut 05155 REAL, DIMENSION(2) :: ta 05156 INTEGER :: i 05157 INTEGER :: iact 05158 INTEGER :: iagain 05159 INTEGER :: idx 05160 INTEGER :: info 05161 INTEGER :: itgbi 05162 INTEGER :: ivgbi 05163 INTEGER :: jcalcm 05164 INTEGER :: k 05165 INTEGER :: labelg 05166 INTEGER :: litera 05167 INTEGER :: lrej 05168 INTEGER :: lun 05169 INTEGER :: lunp 05170 INTEGER :: minf 05171 INTEGER :: mrati 05172 INTEGER :: nan 05173 INTEGER :: nfaci 05174 INTEGER :: nloop 05175 INTEGER :: nrati 05176 INTEGER :: nrej 05177 INTEGER :: nsol 05178 INTEGER :: inone 05179 05180 DOUBLE PRECISION :: stp 05181 DOUBLE PRECISION :: dratio 05182 DOUBLE PRECISION :: dfacin 05183 DOUBLE PRECISION :: djrat 05184 DOUBLE PRECISION :: dwmean 05185 DOUBLE PRECISION :: db 05186 DOUBLE PRECISION :: db1 05187 DOUBLE PRECISION :: db2 05188 DOUBLE PRECISION :: dbdot 05189 LOGICAL :: warner 05190 SAVE 05191 ! ... 05192 05193 ! Printout of algorithm for solution and important parameters ------ 05194 05195 lun=lunlog ! log file 05196 IF(lunlog == 0) lunlog=6 05197 05198 DO lunp=6,lunlog,lunlog-6 05199 WRITE(lunp,*) ' ' 05200 WRITE(lunp,*) 'Solution algorithm: ' 05201 WRITE(lunp,121) '=================================================== ' 05202 05203 IF(metsol == 1) THEN 05204 WRITE(lunp,121) 'solution method:','matrix inversion' 05205 ELSE IF(metsol == 2) THEN 05206 WRITE(lunp,121) 'solution method:','diagonalization' 05207 ELSE 05208 IF(metsol == 3) THEN 05209 WRITE(lunp,121) 'solution method:', 'minres (Paige/Saunders)' 05210 ELSE IF(metsol == 4) THEN 05211 WRITE(lunp,121) 'solution method:', & 05212 'gmres (generalized minimzation of residuals)' 05213 END IF 05214 END IF 05215 WRITE(lunp,123) 'convergence limit at Delta F=',dflim 05216 WRITE(lunp,122) 'maximum number of iterations=',mitera 05217 matrit=MIN(matrit,mitera) 05218 IF(matrit > 1) THEN 05219 WRITE(lunp,122) 'matrix recalculation up to ',matrit, '. iteration' 05220 END IF 05221 IF(metsol >= 3) THEN 05222 IF(matsto == 1) THEN 05223 WRITE(lunp,121) 'matrix storage:','full' 05224 ELSE IF(matsto == 2) THEN 05225 WRITE(lunp,121) 'matrix storage:','sparse' 05226 END IF 05227 WRITE(lunp,122) 'pre-con band-width parameter=',mbandw 05228 IF(mbandw == 0) THEN 05229 WRITE(lunp,121) 'pre-conditioning:','default' 05230 ELSE IF(mbandw < 0) THEN 05231 WRITE(lunp,121) 'pre-conditioning:','none!' 05232 ELSE IF(mbandw > 0) THEN 05233 WRITE(lunp,121) 'pre-conditioning=','band-matrix' 05234 END IF 05235 END IF 05236 IF(regpre == 0.0D0.AND.npresg == 0) THEN 05237 WRITE(lunp,121) 'using pre-sigmas:','no' 05238 ELSE 05239 ! FIXME: NPRESG contains parameters that failed the 'entries' cut... 05240 WRITE(lunp,124) 'pre-sigmas defined for', & 05241 FLOAT(100*npresg)/FLOAT(nvgb),' % of variable parameters' 05242 WRITE(lunp,123) 'default pre-sigma=',regpre 05243 END IF 05244 IF(nregul == 0) THEN 05245 WRITE(lunp,121) 'regularization:','no' 05246 ELSE 05247 WRITE(lunp,121) 'regularization:','yes' 05248 WRITE(lunp,123) 'regularization factor=',regula 05249 END IF 05250 05251 IF(chicut /= 0.0) THEN 05252 WRITE(lunp,121) 'Chi square cut equiv 3 st.dev applied' 05253 WRITE(lunp,123) '... in first iteration with factor',chicut 05254 WRITE(lunp,123) '... in second iteration with factor',chirem 05255 WRITE(lunp,121) ' (reduced by sqrt in next iterations)' 05256 END IF 05257 05258 IF(lhuber /= 0) THEN 05259 WRITE(lunp,122) 'Down-weighting of outliers in', lhuber,' iterations' 05260 WRITE(lunp,123) 'Cut on downweight fraction',dwcut 05261 END IF 05262 05263 05264 121 FORMAT(1X,a40,3X,a) 05265 122 FORMAT(1X,a40,2X,i2,a) 05266 123 FORMAT(1X,a40,2X,e9.2) 05267 124 FORMAT(1X,a40,3X,f5.1,a) 05268 END DO 05269 05270 ! initialization of iterations ------------------------------------- 05271 05272 iitera=0 05273 nloop =0 05274 nsol =0 ! counter for solutions 05275 info =0 05276 lsinfo=0 05277 stp =0.0D0 05278 stepl =SNGL(stp) 05279 concut=1.0E-12 ! initial constraint accuracy 05280 concu2=1.0E-06 ! constraint accuracy 05281 icalcm=1 ! require matrix calculation 05282 iterat=0 ! iteration counter 05283 iterat=-1 05284 litera=-2 05285 nrej=0 ! reset number of rejects 05286 IF(metsol == 1) THEN 05287 wolfc2=0.5 ! not accurate 05288 minf=1 05289 ELSE IF(metsol == 2) THEN 05290 wolfc2=0.5 ! not acurate 05291 minf=2 05292 ELSE IF(metsol == 3) THEN 05293 wolfc2=0.1 ! accurate 05294 minf=3 05295 ELSE IF(metsol == 4) THEN 05296 wolfc2=0.1 ! accurate 05297 minf=3 05298 END IF 05299 05300 ! check initial feasibility of constraint equations ---------------- 05301 05302 WRITE(*,*) ' ' 05303 IF(nofeas == 0) THEN ! make parameter feasible 05304 WRITE(lunlog,*) 'Checking feasibility of parameters:' 05305 WRITE(*,*) 'Checking feasibility of parameters:' 05306 CALL feasib(concut,iact) ! check feasibility 05307 IF(iact /= 0) THEN ! done ... 05308 WRITE(*,102) concut 05309 WRITE(*,*) ' parameters are made feasible' 05310 WRITE(lunlog,102) concut 05311 WRITE(lunlog,*) ' parameters are made feasible' 05312 ELSE ! ... was OK 05313 WRITE(*,*) ' parameters are feasible (i.e. satisfy constraints)' 05314 WRITE(lunlog,*) ' parameters are feasible (i.e. satisfy constraints)' 05315 END IF 05316 concut=concu2 ! cut for constraint check 05317 END IF 05318 iact=1 ! set flag for new data loop 05319 nofeas=0 ! set check-feasibility flag 05320 05321 WRITE(*,*) ' ' 05322 WRITE(*,*)'Reading files and accumulating vectors/matrices ...' 05323 WRITE(*,*) ' ' 05324 05325 CALL etime(ta,rstart) 05326 iterat=-1 05327 litera= 0 05328 jcalcm=-1 05329 iagain= 0 05330 05331 icalcm=1 05332 05333 ! Block 1: data loop with vector (and matrix) calculation ---------- 05334 05335 DO 05336 IF(iterat >= 0) THEN 05337 lcalcm=jcalcm+3 ! mode (1..4) of last loop 05338 IF(jcalcm+1 /= 0) THEN 05339 IF(iterat == 0) THEN 05340 CALL ploopa(6) ! header 05341 CALL ploopb(6) 05342 CALL ploopa(lunlog) ! iteration line 05343 CALL ploopb(lunlog) 05344 iterat=1 05345 CALL gmpxyd(1,FLOAT(nloopn),REAL(fvalue),0.5,0.) ! fcn-value graph (no Delta) 05346 ELSE 05347 IF(iterat /= litera) THEN 05348 CALL ploopb(6) 05349 ! CALL PLOOPA(LUNLOG) 05350 CALL ploopb(lunlog) 05351 litera=iterat 05352 CALL gmpxyd(1,FLOAT(nloopn),REAL(fvalue),0.5,delfun) ! fcn-value (with expected) 05353 IF(metsol == 3) THEN ! extend to 4, i.e. GMRES? 05354 CALL gmpxy(2,FLOAT(iterat),FLOAT(iitera)) ! MINRES iterations 05355 END IF 05356 ELSE 05357 CALL ploopc(6) ! sub-iteration line 05358 CALL ploopc(lunlog) 05359 CALL gmpxyd(1,FLOAT(nloopn),REAL(fvalue),0.5,0.) ! fcn-value graph (no Delta) 05360 END IF 05361 END IF 05362 ELSE 05363 CALL ploopd(6) ! solution line 05364 CALL ploopd(lunlog) 05365 END IF 05366 CALL etime(ta,rstart) 05367 ! CHK 05368 IF (IABS(jcalcm) <= 1) THEN 05369 idx=jcalcm+4 05370 times(idx )=(times(idx )*times(idx+3)+deltim) /(times(idx+3)+1.0) 05371 times(idx+3)= times(idx+3)+1.0 05372 END IF 05373 END IF 05374 jcalcm=icalcm 05375 05376 IF(icalcm >= 0) THEN ! ICALCM = +1 & 0 05377 CALL loopn ! data loop 05378 CALL addcst ! constraints 05379 lrej=nrej 05380 nrej=nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) ! total number of rejects 05381 IF(3*nrej > nrecal) THEN 05382 WRITE(*,*) ' ' 05383 WRITE(*,*) 'Data rejected in previous loop: ' 05384 WRITE(*,*) ' ', & 05385 nrejec(0), ' (rank deficit/NaN) ',nrejec(1),' (Ndf=0) ', & 05386 nrejec(2), ' (huge) ',nrejec(3),' (large)' 05387 WRITE(*,*) 'Too many rejects (>33.3%) - stop' 05388 STOP 05389 END IF 05390 END IF 05391 ! Block 2: new iteration with calculation of solution -------------- 05392 05393 IF(ABS(icalcm) == 1) THEN ! ICALCM = +1 & -1 05394 DO i=1,nagb 05395 globalCorrections(i)=globalVector(i) ! copy rhs 05396 END DO 05397 DO i=1,nvgb 05398 itgbi=globalParVarToTotal(i) 05399 workspaceLinesearch(i)=globalParameter(itgbi) ! copy X for line search 05400 END DO 05401 iterat=iterat+1 ! increase iteration count 05402 IF(metsol == 1) THEN 05403 CALL minver ! inversion 05404 ELSE IF(metsol == 2) THEN 05405 CALL mdiags ! diagonalization 05406 ELSE IF(metsol == 3) THEN 05407 CALL mminrs ! MINRES 05408 ELSE IF(metsol == 4) THEN 05409 WRITE(*,*) '... reserved for GMRES (not yet!)' 05410 CALL mminrs ! GMRES not yet 05411 END IF 05412 05413 ! check feasibility and evtl. make step vector feasible 05414 05415 DO i=1,nvgb 05416 itgbi=globalParVarToTotal(i) 05417 globalParCopy(itgbi)=globalParameter(itgbi) ! save 05418 globalParameter(itgbi)=globalParameter(itgbi)+globalCorrections(i) ! update 05419 END DO 05420 CALL feasib(concut,iact) ! improve constraints 05421 concut=concu2 ! new cut for constraint check 05422 DO i=1,nvgb 05423 itgbi=globalParVarToTotal(i) 05424 globalCorrections(i)=globalParameter(itgbi)-globalParCopy(itgbi) ! feasible stp 05425 globalParameter(itgbi)=globalParCopy(itgbi) ! restore 05426 END DO 05427 05428 ! initialize line search based on slopes and prepare next 05429 05430 CALL ptldef(wolfc2, 10.0, minf,10) 05431 IF(metsol == 1) THEN 05432 wolfc2=0.5 ! not accurate 05433 minf=3 05434 ELSE IF(metsol == 2) THEN 05435 wolfc2=0.5 ! not acurate 05436 minf=3 05437 ELSE IF(metsol == 3) THEN 05438 wolfc2=0.1 ! accurate 05439 minf=4 05440 ELSE IF(metsol == 4) THEN 05441 wolfc2=0.1 ! accurate 05442 minf=4 05443 END IF 05444 db=dbdot(nvgb,globalCorrections,globalVector) 05445 db1=dbdot(nvgb,globalCorrections,globalCorrections) 05446 db2=dbdot(nvgb,globalVector,globalVector) 05447 delfun=SNGL(db) 05448 angras=SNGL(db/DSQRT(db1*db2)) 05449 05450 IF(db <= 0.0D0) THEN 05451 WRITE(*,*) 'Function not decreasing:',db 05452 IF(db <= -1.0D-3) THEN ! 100311, VB/CK: allow some margin for numerics 05453 iagain=iagain+1 05454 IF (iagain <= 1) THEN 05455 WRITE(*,*) '... again matrix calculation' 05456 icalcm=1 05457 CYCLE 05458 ELSE 05459 WRITE(*,*) '... aborting iterations' 05460 GO TO 90 05461 END IF 05462 ELSE 05463 WRITE(*,*) '... stopping iterations' 05464 iagain=0 05465 GO TO 90 05466 END IF 05467 ELSE 05468 iagain=0 05469 END IF 05470 icalcm=0 ! switch 05471 ENDIF 05472 ! Block 3: line searching ------------------------------------------ 05473 05474 IF(icalcm+2 == 0) EXIT 05475 CALL ptline(nagb,workspaceLinesearch, & ! current parameter values 05476 flines, & ! chi^2 function value 05477 globalVector, & ! gradient 05478 globalCorrections, & ! step vector stp 05479 stp, & ! returned step factor 05480 info) ! returned information 05481 ! WRITE(*,*) 'PTLINE returns INFO, STP=',INFO, STP 05482 lsinfo=info 05483 05484 stepl=SNGL(stp) 05485 nan=0 05486 DO i=1,nvgb 05487 itgbi=globalParVarToTotal(i) 05488 IF ((.NOT.(workspaceLinesearch(i) <= 0.0D0)).AND. & 05489 (.NOT.(workspaceLinesearch(i) > 0.0D0))) nan=nan+1 05490 globalParameter(itgbi)=workspaceLinesearch(i) ! current parameter values 05491 END DO 05492 05493 IF (nan > 0) THEN 05494 WRITE(*,*) 'Result vector containes ', nan,' NaNs - stop' 05495 STOP 05496 END IF 05497 05498 ! subito exit, if required ----------------------------------------- 05499 05500 IF(isubit /= 0) THEN ! subito 05501 WRITE(*,*) 'Subito! Exit after first step.' 05502 GO TO 90 05503 END IF 05504 05505 IF(info == 0) THEN 05506 WRITE(*,*) 'INFO=0 should not happen (line search input err)' 05507 IF (iagain <= 0) THEN 05508 icalcm=1 05509 CYCLE 05510 ENDIF 05511 END IF 05512 IF(info < 0) CYCLE 05513 ! Block 4: line search convergence --------------------------------- 05514 05515 CALL ptlprt(lunlog) 05516 CALL feasib(concut,iact) ! check constraints 05517 IF(iact /= 0.OR.chicut > 1.0) THEN 05518 icalcm=-1 05519 IF(iterat < matrit) icalcm=+1 05520 CYCLE ! iterate 05521 END IF 05522 IF(delfun <= dflim) GO TO 90 ! convergence 05523 IF(iterat >= mitera) GO TO 90 ! ending 05524 icalcm=-1 05525 IF(iterat < matrit) icalcm=+1 05526 CYCLE ! next iteration 05527 05528 ! Block 5: iteration ending ---------------------------------------- 05529 05530 90 icalcm=-2 05531 END DO 05532 IF(nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) /= 0) THEN 05533 WRITE(*,*) ' ' 05534 WRITE(*,*) 'Data rejected in last loop: ' 05535 WRITE(*,*) ' ', & 05536 nrejec(0), ' (rank deficit/NaN) ',nrejec(1),' (Ndf=0) ', & 05537 nrejec(2), ' (huge) ',nrejec(3),' (large)' 05538 END IF 05539 05540 dwmean=sumndf/dfloat(ndfsum) 05541 dratio=fvalue/dwmean/dfloat(ndfsum-nagb) 05542 catio=SNGL(dratio) 05543 IF(lhuber /= 0) THEN 05544 IF (lhuber <= 3) THEN 05545 catio=catio/0.9326 ! correction Huber downweighting 05546 ELSE 05547 catio=catio/0.8228 ! correction Cauchy downweighting 05548 END IF 05549 END IF 05550 mrati=nint(100.0*catio) 05551 05552 DO lunp=6,lunlog,lunlog-6 05553 WRITE(lunp,*) ' ' 05554 IF (nfilw <= 0) THEN 05555 WRITE(lunp,*) 'Sum(Chi^2)/Sum(Ndf) =',fvalue 05556 WRITE(lunp,*) ' / (',ndfsum,'-',nagb,')' 05557 WRITE(lunp,*) ' =',dratio 05558 ELSE 05559 WRITE(lunp,*) 'Sum(W*Chi^2)/Sum(Ndf)/<W> =',fvalue 05560 WRITE(lunp,*) ' / (',ndfsum,'-', nagb,')' 05561 WRITE(lunp,*) ' /',dwmean 05562 WRITE(lunp,*) ' =',dratio 05563 END IF 05564 WRITE(lunp,*) ' ' 05565 IF(lhuber /= 0) WRITE(lunp,*) & 05566 ' with correction for down-weighting ',catio 05567 END DO 05568 nrej=nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) ! total number of rejects 05569 05570 ! ... the end with exit code ??????????????????????????????????????? 05571 05572 ! WRITE(*,199) ! write exit code 05573 ! + '-----------------------------------------------------------' 05574 ! IF(ITEXIT.EQ.0) WRITE(*,199) 05575 ! + 'Exit code = 0: Convergence reached' 05576 ! IF(ITEXIT.EQ.1) WRITE(*,199) 05577 ! + 'Exit code = 1: No improvement in last iteration' 05578 ! IF(ITEXIT.EQ.2) WRITE(*,199) 05579 ! + 'Exit code = 2: Maximum number of iterations reached' 05580 ! IF(ITEXIT.EQ.3) WRITE(*,199) 05581 ! + 'Exit code = 3: Failure' 05582 ! WRITE(*,199) 05583 ! + '-----------------------------------------------------------' 05584 ! WRITE(*,199) ' ' 05585 05586 05587 nrati=nint(10000.0*FLOAT(nrej)/FLOAT(nrecal)) 05588 djrat=0.01D0*dfloat(nrati) 05589 nfaci=nint(100.0*SQRT(catio)) 05590 05591 dratio=0.01D0*dfloat(mrati) 05592 dfacin=0.01D0*dfloat(nfaci) 05593 05594 warner=.FALSE. 05595 IF(mrati < 90.OR.mrati > 110) warner=.TRUE. 05596 IF(nrati > 1) warner=.TRUE. 05597 IF(nmiss1 /= 0) warner=.TRUE. 05598 IF(iagain /= 0) warner=.TRUE. 05599 05600 IF(warner) THEN 05601 WRITE(*,199) ' ' 05602 WRITE(*,199) ' ' 05603 WRITE(*,199) 'WarningWarningWarningWarningWarningWarningWarningWarningWar' 05604 WRITE(*,199) 'arningWarningWarningWarningWarningWarningWarningWarningWarn' 05605 WRITE(*,199) 'rningWarningWarningWarningWarningWarningWarningWarningWarni' 05606 WRITE(*,199) 'ningWarningWarningWarningWarningWarningWarningWarningWarnin' 05607 WRITE(*,199) 'ingWarningWarningWarningWarningWarningWarningWarningWarning' 05608 WRITE(*,199) 'ngWarningWarningWarningWarningWarningWarningWarningWarningW' 05609 WRITE(*,199) 'gWarningWarningWarningWarningWarningWarningWarningWarningWa' 05610 05611 IF(mrati < 90.OR.mrati > 110) THEN 05612 WRITE(*,199) ' ' 05613 WRITE(*,*) ' Chi^2/Ndf = ',dratio, ' (should be close to 1)' 05614 WRITE(*,*) ' => multiply all input standard ', & 05615 'deviations by factor',dfacin 05616 END IF 05617 05618 IF(nrati > 1) THEN 05619 WRITE(*,199) ' ' 05620 WRITE(*,*) ' Fraction of rejects =',djrat,' %', & 05621 ' (should be far below 1 %)' 05622 WRITE(*,*) ' => please provide correct mille data' 05623 END IF 05624 05625 IF(iagain /= 0) THEN 05626 WRITE(*,199) ' ' 05627 WRITE(*,*) ' Matrix not positiv definite '// & 05628 '(function not decreasing)' 05629 WRITE(*,*) ' => please provide correct mille data' 05630 END IF 05631 05632 IF(nmiss1 /= 0) THEN 05633 WRITE(*,199) ' ' 05634 WRITE(*,*) ' Rank defect =',nmiss1, & 05635 ' for constraint equations, should be 0' 05636 WRITE(*,*) ' => please correct constraint definition' 05637 END IF 05638 05639 WRITE(*,199) ' ' 05640 WRITE(*,199) 'WarningWarningWarningWarningWarningWarningWarningWarningWar' 05641 WRITE(*,199) 'arningWarningWarningWarningWarningWarningWarningWarningWarn' 05642 WRITE(*,199) 'rningWarningWarningWarningWarningWarningWarningWarningWarni' 05643 WRITE(*,199) 'ningWarningWarningWarningWarningWarningWarningWarningWarnin' 05644 WRITE(*,199) 'ingWarningWarningWarningWarningWarningWarningWarningWarning' 05645 WRITE(*,199) 'ngWarningWarningWarningWarningWarningWarningWarningWarningW' 05646 WRITE(*,199) 'gWarningWarningWarningWarningWarningWarningWarningWarningWa' 05647 WRITE(*,199) ' ' 05648 05649 ENDIF 05650 05651 CALL mend ! modul ending 05652 05653 ! ------------------------------------------------------------------ 05654 05655 IF(metsol == 1) THEN 05656 05657 ELSE IF(metsol == 2) THEN 05658 CALL zdiags 05659 ELSE IF(metsol == 3) THEN 05660 ! errors and correlations from MINRES 05661 DO k=1,mnrsel 05662 labelg=lbmnrs(k) 05663 IF(labelg == 0) CYCLE 05664 itgbi=inone(labelg) 05665 ivgbi=0 05666 IF(itgbi /= 0) ivgbi=globalParLabelIndex(2,itgbi) 05667 IF(ivgbi < 0) ivgbi=0 05668 IF(ivgbi == 0) CYCLE 05669 ! determine error and global correlation for parameter IVGBI 05670 CALL solglo(ivgbi) 05671 END DO 05672 05673 ELSE IF(metsol == 4) THEN 05674 05675 END IF 05676 05677 CALL prtglo ! print result 05678 05679 102 FORMAT(' Call FEASIB with cut=',g10.3) 05680 ! 103 FORMAT(1X,A,G12.4) 05681 199 FORMAT(7X,a) 05682 END SUBROUTINE xloopn ! standard solution 05683 05697 05698 SUBROUTINE filetc 05699 USE mpmod 05700 USE mpdalc 05701 05702 IMPLICIT NONE 05703 INTEGER :: i 05704 INTEGER :: ia 05705 INTEGER :: iargc 05706 INTEGER :: ib 05707 INTEGER :: ie 05708 INTEGER :: ierrf 05709 INTEGER :: ieq 05710 INTEGER :: ioff 05711 INTEGER :: iopt 05712 INTEGER :: ios 05713 INTEGER :: iosum 05714 INTEGER :: it 05715 INTEGER :: k 05716 INTEGER :: mat 05717 INTEGER :: nab 05718 INTEGER :: nline 05719 INTEGER :: npat 05720 INTEGER :: ntext 05721 INTEGER :: nu 05722 INTEGER :: nuf 05723 INTEGER :: nums 05724 INTEGER :: nufile 05725 INTEGER :: lenfileInfo 05726 INTEGER :: lenFileNames 05727 INTEGER :: matint 05728 INTEGER, DIMENSION(:,:), ALLOCATABLE :: vecfileInfo 05729 INTEGER, DIMENSION(:,:), ALLOCATABLE :: tempArray 05730 INTEGER(KIND=large) :: rows 05731 INTEGER(KIND=large) :: cols 05732 INTEGER(KIND=large) :: newcols 05733 INTEGER(KIND=large) :: length 05734 05735 CHARACTER (LEN=1024) :: text 05736 CHARACTER (LEN=1024) :: fname 05737 CHARACTER (LEN=14) :: bite(3) 05738 CHARACTER (LEN=32) :: keystx 05739 DOUBLE PRECISION :: dnum(100) 05740 SAVE 05741 DATA bite/'C_binary','text ','Fortran_binary'/ 05742 ! ... 05743 CALL mstart('FILETC/X') 05744 05745 nuf=1 ! C binary is default 05746 DO i=1,8 05747 times(i)=0.0 05748 END DO 05749 05750 ! read command line options ---------------------------------------- 05751 05752 filnam=' ' ! print command line options and find steering file 05753 DO i=1,iargc() 05754 IF(i == 1) THEN 05755 WRITE(*,*) ' ' 05756 WRITE(*,*) 'Command line options: ' 05757 WRITE(*,*) '--------------------- ' 05758 END IF 05759 CALL getarg(i,text) ! get I.th text from command line 05760 CALL rltext(text,ia,ib,nab) ! return indices for non-blank area 05761 WRITE(*,101) i,text(1:nab) ! echo print 05762 IF(text(ia:ia) /= '-') THEN 05763 nu=nufile(text(ia:ib)) ! inquire on file existence 05764 IF(nu == 2) THEN ! existing text file 05765 IF(filnam /= ' ') THEN 05766 WRITE(*,*) 'Second text file in command line - stop' 05767 STOP 05768 ELSE 05769 filnam=text 05770 END IF 05771 ELSE 05772 WRITE(*,*) 'Open error for file:',text(ia:ib),' - stop' 05773 STOP 05774 END IF 05775 ELSE 05776 IF(INDEX(text(ia:ib),'b') /= 0) THEN 05777 mdebug=3 ! debug flag 05778 WRITE(*,*) 'Debugging requested' 05779 END IF 05780 it=INDEX(text(ia:ib),'t') 05781 IF(it /= 0) THEN 05782 ictest=1 ! internal test files 05783 ieq=INDEX(text(ia+it:ib),'=')+it 05784 IF (it /= ieq) THEN 05785 IF (INDEX(text(ia+ieq:ib),'SL0' ) /= 0) ictest=2 05786 IF (INDEX(text(ia+ieq:ib),'SLE' ) /= 0) ictest=3 05787 IF (INDEX(text(ia+ieq:ib),'BP' ) /= 0) ictest=4 05788 IF (INDEX(text(ia+ieq:ib),'BRLF') /= 0) ictest=5 05789 IF (INDEX(text(ia+ieq:ib),'BRLC') /= 0) ictest=6 05790 END IF 05791 END IF 05792 IF(INDEX(text(ia:ib),'s') /= 0) isubit=1 ! like "subito" 05793 IF(INDEX(text(ia:ib),'f') /= 0) iforce=1 ! like "force" 05794 END IF 05795 IF(i == iargc()) WRITE(*,*) '--------------------- ' 05796 END DO 05797 05798 05799 ! create test files for option -t ---------------------------------- 05800 05801 IF(ictest >= 1) THEN 05802 WRITE(*,*) ' ' 05803 IF (ictest == 1) THEN 05804 CALL mptest ! 'wire chamber' 05805 ELSE 05806 CALL mptst2(ictest-2) ! 'silicon tracker' 05807 END IF 05808 IF(filnam == ' ') filnam='mp2str.txt' 05809 WRITE(*,*) ' ' 05810 END IF 05811 05812 ! check default steering file with file-name "steerfile" ----------- 05813 05814 IF(filnam == ' ') THEN ! check default steering file 05815 text='steerfile' 05816 CALL rltext(text,ia,ib,nab) ! return indices for non-blank area 05817 nu=nufile(text(ia:ib)) ! inquire on file existence and type 05818 IF(nu > 0) THEN 05819 filnam=text 05820 ELSE 05821 STOP 'in FILETC: no steering file. .' 05822 END IF 05823 END IF 05824 05825 05826 ! open, read steering file: 05827 ! end 05828 ! fortranfiles 05829 ! cfiles 05830 05831 05832 CALL rltext(filnam,ia,ib,nfnam) ! return indices for non-blank area 05833 WRITE(*,*) ' ' 05834 WRITE(*,*) 'Listing of steering file: ',filnam(1:nfnam) 05835 WRITE(*,*) '-------------------------' 05836 OPEN(10,FILE=filnam(1:nfnam),IOSTAT=ios) 05837 IF(ios /= 0) THEN 05838 WRITE(*,*) 'Open error for steering file - stop' 05839 STOP 05840 END IF 05841 ifile =0 05842 nfiles=0 05843 05844 lenfileInfo=2 05845 lenFileNames=0 05846 rows=6; cols=lenFileInfo 05847 CALL mpalloc(vecfileInfo,rows,cols,'file info from steering') 05848 nline=0 05849 DO 05850 READ(10,102,IOSTAT=ierrf) text ! read steering file 05851 IF (ierrf < 0) EXIT ! eof 05852 CALL rltext(text,ia,ib,nab) ! return indices for non-blank area 05853 nline=nline+1 05854 IF(nline <= 50) THEN ! print up to 50 lines 05855 WRITE(*,101) nline,text(1:nab) 05856 IF(nline == 50) WRITE(*,*) ' ...' 05857 END IF 05858 05859 CALL rltext(text,ia,ib,nab) ! test content 'end' 05860 IF(ib == ia+2) THEN 05861 mat=matint(text(ia:ib),'end',npat,ntext) 05862 IF(mat == 3) THEN 05863 text=' ' 05864 CALL intext(text,nline) 05865 WRITE(*,*) ' end-statement after',nline,' text lines' 05866 EXIT 05867 END IF 05868 END IF 05869 05870 keystx='fortranfiles' 05871 !GF MAT=MATINT(TEXT(IA:IB),KEYSTX,NPAT,NTEXT) 05872 !GF IF(MAT.GE.NTEXT-NTEXT/10) THEN 05873 IF (text(ia:ib) == keystx) THEN ! exact comparison by GF 05874 nuf=3 05875 ! WRITE(*,*) 'Fortran files' 05876 CYCLE 05877 END IF 05878 05879 keystx='Cfiles' 05880 !GF MAT=MATINT(TEXT(IA:IB),KEYSTX,NPAT,NTEXT) 05881 !GF IF(MAT.GE.NTEXT-NTEXT/10) THEN 05882 IF (text(ia:ib) == keystx) THEN ! exact comparison by GF 05883 nuf=1 05884 ! WRITE(*,*) 'Cfiles' 05885 CYCLE 05886 END IF 05887 05888 ! file names 05889 ! check for file options (' -- ') 05890 ie=ib 05891 iopt=INDEX(text(ia:ib),' -- ') 05892 IF (iopt > 0) ie=iopt-1 05893 05894 IF(nab == 0) CYCLE 05895 nu=nufile(text(ia:ie)) ! inquire on file existence 05896 IF(nu > 0) THEN ! existing file 05897 IF (nfiles == lenFileInfo) THEN ! increase length 05898 CALL mpalloc(tempArray,rows,cols,'temp file info from steering') 05899 tempArray=vecfileInfo 05900 CALL mpdealloc(vecfileInfo) 05901 lenFileInfo=lenFileInfo*2 05902 newcols=lenFileInfo 05903 CALL mpalloc(vecfileInfo,rows,newcols,'file info from steering') 05904 vecfileInfo(:,1:cols)=tempArray(:,1:cols) 05905 CALL mpdealloc(tempArray) 05906 cols=newcols 05907 ENDIF 05908 nfiles=nfiles+1 ! count number of files 05909 IF(nu == 1) nu=nuf ! 05910 lenFileNames=lenFileNames+ie-ia+1 ! total length of file names 05911 vecFileInfo(1,nfiles)=nline ! line number 05912 vecFileInfo(2,nfiles)=nu ! cbinary =1, text =2, fbinary=3 05913 vecFileInfo(3,nfiles)=ia ! file name start 05914 vecFileInfo(4,nfiles)=ie ! file name end 05915 vecFileInfo(5,nfiles)=iopt ! option start 05916 vecFileInfo(6,nfiles)=ib ! option end 05917 ELSE 05918 ! WRITE(*,*) 'Open error for file ',TEXT(IA:IB) 05919 ! STOP 05920 END IF 05921 END DO 05922 REWIND 10 05923 ! read again to fill dynamic arrays with file info 05924 length=nfiles 05925 CALL mpalloc(mfd,length,'file type') 05926 CALL mpalloc(nfd,length,'file line (in steering)') 05927 CALL mpalloc(lfd,length,'file name length') 05928 CALL mpalloc(ofd,length,'file option') 05929 length=lenFileNames 05930 CALL mpalloc(tfd,length,'file name') 05931 nline=0 05932 i=1 05933 ioff=0 05934 DO 05935 READ(10,102,IOSTAT=ierrf) text ! read steering file 05936 IF (ierrf < 0) EXIT ! eof 05937 nline=nline+1 05938 IF (nline == vecFileInfo(1,i)) THEN 05939 nfd(i)=vecFileInfo(1,i) 05940 mfd(i)=vecFileInfo(2,i) 05941 ia=vecFileInfo(3,i)-1 05942 lfd(i)=vecFileInfo(4,i)-ia ! length file name 05943 FORALL (k=1:lfd(i)) tfd(ioff+k)=text(ia+k:ia+k) 05944 ! tfd(i)=text(vecFileInfo(3,i):vecFileInfo(4,i)) ! file name 05945 ioff=ioff+lfd(i) 05946 ofd(i)=1.0 ! option for file 05947 IF (vecFileInfo(5,i) > 0) THEN 05948 CALL ratext(text(vecFileInfo(5,i)+4:vecFileInfo(6,i)),nums,dnum) ! translate text to DP numbers 05949 IF (nums > 0) ofd(i)=SNGL(dnum(1)) 05950 END IF 05951 i=i+1 05952 IF (i > nfiles) EXIT 05953 ENDIF 05954 ENDDO 05955 CALL mpdealloc(vecfileInfo) 05956 REWIND 10 05957 ! additional info for binary files 05958 length=nfiles; rows=2 05959 CALL mpalloc(ifd,length,'integrated record numbers (=offset)') 05960 CALL mpalloc(jfd,length,'number of accepted records') 05961 CALL mpalloc(kfd,rows,length,'number of records in file, file order') 05962 CALL mpalloc(dfd,length,'ndf sum') 05963 CALL mpalloc(xfd,length,'max. record size') 05964 CALL mpalloc(wfd,length,'file weight') 05965 CALL mpalloc(cfd,length,'chi2 sum') 05966 ! 05967 WRITE(*,*) '-------------------------' 05968 WRITE(*,*) ' ' 05969 05970 ! print table of files --------------------------------------------- 05971 05972 IF (mprint > 1) THEN 05973 WRITE(*,*) 'Table of files:' 05974 WRITE(*,*) '---------------' 05975 END IF 05976 WRITE(8,*) ' ' 05977 WRITE(8,*) 'Text and data files:' 05978 ioff=0 05979 DO i=1,nfiles 05980 FORALL (k=1:lfd(i)) fname(k:k)=tfd(ioff+k) 05981 ! fname=tfd(i)(1:lfd(i)) 05982 IF (mprint > -1) WRITE(*,103) i,bite(mfd(i)),fname(1:lfd(i)) 05983 WRITE(8,103) i,bite(mfd(i)),fname(1:lfd(i)) 05984 ioff=ioff+lfd(i) 05985 END DO 05986 IF (mprint > 1) THEN 05987 WRITE(*,*) '---------------' 05988 WRITE(*,*) ' ' 05989 END IF 05990 05991 ! open the binary (data) files on unit 11, 12, ... 05992 05993 iosum=0 05994 nfilf=0 05995 nfilb=0 05996 nfilw=0 05997 ioff=0 05998 DO i=1,nfiles ! Fortran files 05999 IF(mfd(i) == 3) THEN 06000 FORALL (k=1:lfd(i)) fname(k:k)=tfd(ioff+k) 06001 ! fname=tfd(i)(1:lfd(i)) 06002 WRITE(*,*) 'Opening Fortran file ',10+nfilf+1, ' ',fname(1:lfd(i)) 06003 OPEN(10+nfilf+1,FILE=fname(1:lfd(i)),IOSTAT=ios, FORM='UNFORMATTED') 06004 IF(ios /= 0) THEN 06005 WRITE(*,*) 'Open error for file ',fname(1:lfd(i)) 06006 iosum=iosum+1 06007 ELSE 06008 nfilf=nfilf+1 06009 nfilb=nfilb+1 06010 wfd(nfilb)=ofd(i) 06011 END IF 06012 END IF 06013 ioff=ioff+lfd(i) 06014 END DO 06015 06016 ! open the binary C files 06017 06018 nfilc=-1 06019 ioff=0 06020 DO i=1,nfiles ! Cfiles 06021 IF(mfd(i) == 1) THEN 06022 #ifdef READ_C_FILES 06023 IF(nfilc < 0) CALL initc(nfiles) ! uncommented by GF 06024 IF(nfilc < 0) nfilc=0 06025 FORALL (k=1:lfd(i)) fname(k:k)=tfd(ioff+k) 06026 ! fname=tfd(i)(1:lfd(i)) 06027 WRITE(*,*) 'Opening C file ',nfilc+1, ': ',fname(1:lfd(i)) ! by GF 06028 CALL openc(fname(1:lfd(i)),ios) 06029 IF(ios /= 0) THEN 06030 WRITE(*,*) 'Open error ',ios,' for file ',fname(1:lfd(i)) 06031 iosum=iosum+1 ! typo fixed by GF 06032 ELSE 06033 nfilc=nfilc+1 06034 nfilb=nfilb+1 06035 wfd(nfilb)=ofd(i) 06036 END IF 06037 #else 06038 WRITE(*,*) 'Opening of C-files not supported.' 06039 ! GF add 06040 iosum=iosum+1 06041 ! GF add end 06042 #endif 06043 END IF 06044 ioff=ioff+lfd(i) 06045 END DO 06046 06047 DO k=1,nfilb 06048 kfd(1,k)=0 ! reset record counters 06049 kfd(2,k)=k ! set file number 06050 ifd(k)=0 ! reset integrated record numbers 06051 xfd(k)=0 ! reset max record size 06052 END DO 06053 06054 IF(iosum /= 0) THEN 06055 STOP 'FILETC: open error ' 06056 END IF 06057 IF(nfilb == 0) THEN 06058 STOP 'FILETC: no binary files ' 06059 END IF 06060 WRITE(*,*) nfilb,' binary files opened' ! corrected by GF 06061 101 FORMAT(i3,2X,a) 06062 102 FORMAT(a) 06063 103 FORMAT(i3,2X,a14,3X,a) 06064 ! CALL mend 06065 RETURN 06066 END SUBROUTINE filetc 06067 06118 06119 SUBROUTINE filetx ! --------------------------------------------------- 06120 USE mpmod 06121 06122 IMPLICIT NONE 06123 INTEGER :: i 06124 INTEGER :: ia 06125 INTEGER :: ib 06126 INTEGER :: ierrf 06127 INTEGER :: ioff 06128 INTEGER :: ios 06129 INTEGER :: iosum 06130 INTEGER :: k 06131 INTEGER :: mat 06132 INTEGER :: nab 06133 INTEGER :: nfiln 06134 INTEGER :: nline 06135 INTEGER :: nlinmx 06136 INTEGER :: npat 06137 INTEGER :: ntext 06138 INTEGER :: matint 06139 06140 ! CALL MSTART('FILETX') 06141 06142 CHARACTER (LEN=1024) :: text 06143 CHARACTER (LEN=1024) :: fname 06144 06145 WRITE(*,*) ' ' 06146 WRITE(*,*) 'Processing text files ...' 06147 WRITE(*,*) ' ' 06148 06149 iosum=0 06150 ioff=0 06151 DO i=0,nfiles 06152 IF(i == 0) THEN 06153 WRITE(*,*) 'File ',filnam(1:nfnam) 06154 nlinmx=100 06155 ELSE 06156 nlinmx=10 06157 ia=ioff 06158 ioff=ioff+lfd(i) 06159 IF(mfd(i) /= 2) CYCLE ! exclude binary files 06160 FORALL (k=1:lfd(i)) fname(k:k)=tfd(ia+k) 06161 WRITE(*,*) 'File ',fname(1:lfd(i)) 06162 IF (mprint > 1) WRITE(*,*) ' ' 06163 OPEN(10,FILE=fname(1:lfd(i)),IOSTAT=ios,FORM='FORMATTED') 06164 IF(ios /= 0) THEN 06165 WRITE(*,*) 'Open error for file ',fname(1:lfd(i)) 06166 iosum=iosum+1 06167 CYCLE 06168 END IF 06169 END IF 06170 06171 nline=0 06172 nfiln=1 06173 ! read text file 06174 DO 06175 READ(10,102,IOSTAT=ierrf) text 06176 IF (ierrf < 0) THEN 06177 text=' ' 06178 CALL intext(text,nline) 06179 WRITE(*,*) ' end-of-file after',nline,' text lines' 06180 EXIT ! eof 06181 ENDIF 06182 nline=nline+1 06183 IF(nline <= nlinmx.AND.mprint > 1) THEN ! print first 10 lines of every text fiLE 06184 CALL rltext(text,ia,ib,nab) 06185 nab=MAX(1,nab) 06186 WRITE(*,101) nline,text(1:nab) 06187 IF(nline == nlinmx) WRITE(*,*) ' ...' 06188 END IF 06189 06190 CALL rltext(text,ia,ib,nab) ! test content 'end' 06191 IF(ib == ia+2) THEN 06192 mat=matint(text(ia:ib),'end',npat,ntext) 06193 IF(mat == 3) THEN 06194 text=' ' 06195 CALL intext(text,nline) 06196 WRITE(*,*) ' end-statement after',nline,' text lines' 06197 EXIT 06198 END IF 06199 END IF 06200 06201 IF(i == 0) THEN ! first text file - exclude lines with file names 06202 IF(nfiln <= nfiles.AND.nline == nfd(nfiln)) THEN 06203 nfiln=nfiln+1 06204 text=' ' 06205 ! WRITE(*,*) 'line is excluded ',TEXT(1:10) 06206 END IF 06207 END IF 06208 ! WRITE(*,*) TEXT(1:40),' < interprete text' 06209 CALL intext(text,nline) ! interprete text 06210 END DO 06211 WRITE(*,*) ' ' 06212 REWIND 10 06213 CLOSE(UNIT=10) 06214 END DO 06215 06216 IF(iosum /= 0) THEN 06217 STOP 'FILETX: open error(s) in text files ' 06218 END IF 06219 06220 WRITE(*,*) '... end of text file processing.' 06221 WRITE(*,*) ' ' 06222 06223 IF(lunkno /= 0) THEN 06224 WRITE(*,*) ' ' 06225 WRITE(*,*) lunkno,' unknown keywords in steering files, ', & 06226 'or file non-existing,' 06227 WRITE(*,*) ' see above!' 06228 WRITE(*,*) '------------> stop' 06229 WRITE(*,*) ' ' 06230 STOP 06231 END IF 06232 06233 ! check methods 06234 06235 IF(metsol == 0) THEN ! if undefined 06236 IF(matsto == 0) THEN ! if undefined 06237 ! METSOL=1 ! default is matrix inversion 06238 ! MATSTO=1 ! default is symmetric matrix 06239 ELSE IF(matsto == 1) THEN ! if symmetric 06240 metsol=3 ! MINRES 06241 ELSE IF(matsto == 2) THEN ! if sparse 06242 metsol=3 ! MINRES 06243 END IF 06244 ELSE IF(metsol == 1) THEN ! if inversion 06245 matsto=1 ! 06246 ELSE IF(metsol == 2) THEN ! if diagonalization 06247 matsto=1 06248 ELSE IF(metsol == 3) THEN ! if MINRES 06249 ! MATSTO=2 or 1 06250 ELSE IF(metsol == 4) THEN ! if GMRES 06251 ! MATSTO=2 or 1 06252 ELSE 06253 WRITE(*,*) 'MINRES forced with sparse matrix!' 06254 WRITE(*,*) ' ' 06255 WRITE(*,*) 'MINRES forced with sparse matrix!' 06256 WRITE(*,*) ' ' 06257 WRITE(*,*) 'MINRES forced with sparse matrix!' 06258 metsol=3 ! forced 06259 matsto=2 ! forced 06260 END IF 06261 IF(matsto > 2) THEN 06262 WRITE(*,*) 'MINRES forced with sparse matrix!' 06263 WRITE(*,*) ' ' 06264 WRITE(*,*) 'MINRES forced with sparse matrix!' 06265 WRITE(*,*) ' ' 06266 WRITE(*,*) 'MINRES forced with sparse matrix!' 06267 metsol=3 ! forced 06268 matsto=2 ! forced 06269 END IF 06270 06271 ! print information about methods and matrix storage modes 06272 06273 WRITE(*,*) ' ' 06274 WRITE(*,*) 'Solution method and matrix-storage mode:' 06275 IF(metsol == 1) THEN 06276 WRITE(*,*) ' METSOL = 1: matrix inversion' 06277 ELSE IF(metsol == 2) THEN 06278 WRITE(*,*) ' METSOL = 2: diagonalization' 06279 ELSE IF(metsol == 3) THEN 06280 WRITE(*,*) ' METSOL = 3: MINRES' 06281 ELSE IF(metsol == 4) THEN 06282 WRITE(*,*) ' METSOL = 4: GMRES (-> MINRES)' 06283 06284 END IF 06285 06286 WRITE(*,*) ' with',mitera,' iterations' 06287 06288 IF(matsto == 1) THEN 06289 WRITE(*,*) ' MATSTO = 1: symmetric matrix, ', '(n*n+n)/2 elements' 06290 ELSE IF(matsto == 2) THEN 06291 WRITE(*,*) ' MATSTO = 2: sparse matrix' 06292 END IF 06293 IF(mbandw /= 0) THEN 06294 WRITE(*,*) ' and band matrix, width',mbandw 06295 END IF 06296 06297 06298 06299 IF(chicut /= 0.0) THEN 06300 WRITE(*,*) 'Chi square cut equiv 3 st.dev applied ...' 06301 WRITE(*,*) ' in first iteration with factor',chicut 06302 WRITE(*,*) ' in second iteration with factor',chirem 06303 WRITE(*,*) ' (reduced by sqrt in next iterations)' 06304 END IF 06305 06306 IF(lhuber /= 0) THEN 06307 WRITE(*,*) ' Down-weighting of outliers in', lhuber,' iterations' 06308 WRITE(*,*) ' Cut on downweight fraction',dwcut 06309 END IF 06310 06311 CALL mend 06312 06313 101 FORMAT(i3,2X,a) 06314 102 FORMAT(a) 06315 END SUBROUTINE filetx 06316 06326 06327 INTEGER FUNCTION nufile(fname) 06328 06329 IMPLICIT NONE 06330 INTEGER :: ios 06331 INTEGER :: l1 06332 INTEGER :: ll 06333 INTEGER :: nm 06334 INTEGER :: npat 06335 INTEGER :: ntext 06336 INTEGER :: nuprae 06337 INTEGER :: matint 06338 06339 CHARACTER (LEN=*), INTENT(INOUT) :: fname 06340 LOGICAL :: ex 06341 SAVE 06342 ! ... 06343 nufile=0 06344 IF(fname(1:5) == 'rfio:') nuprae=1 06345 IF(fname(1:5) == 'dcap:') nuprae=2 06346 IF(fname(1:5) == 'root:') nuprae=3 06347 IF(nuprae == 0) THEN 06348 INQUIRE(FILE=fname,IOSTAT=ios,EXIST=ex) 06349 IF(ios /= 0) nufile=-ABS(ios) 06350 IF(ios /= 0) RETURN 06351 ELSE IF(nuprae == 1) THEN ! rfio: 06352 ll=LEN(fname) 06353 fname=fname(6:ll) 06354 ex=.TRUE. 06355 nufile=1 06356 RETURN 06357 ELSE 06358 ex=.TRUE. ! assume file existence 06359 END IF 06360 IF(ex) THEN 06361 nufile=1 ! binary 06362 ll=LEN(fname) 06363 l1=MAX(1,ll-3) 06364 nm=matint('xt',fname(l1:ll),npat,ntext) 06365 IF(nm == 2) nufile=2 ! text 06366 IF(nm < 2) THEN 06367 nm=matint('tx',fname(l1:ll),npat,ntext) 06368 IF(nm == 2) nufile=2 ! text 06369 END IF 06370 END IF 06371 END FUNCTION nufile 06372 06380 SUBROUTINE intext(text,nline) 06381 USE mpmod 06382 USE mptext 06383 06384 IMPLICIT NONE 06385 INTEGER :: i 06386 INTEGER :: ia 06387 INTEGER :: ib 06388 INTEGER :: icount 06389 INTEGER :: ier 06390 INTEGER :: iomp 06391 INTEGER :: k 06392 INTEGER :: kkey 06393 INTEGER :: label 06394 INTEGER :: lkey 06395 INTEGER :: mat 06396 INTEGER :: miter 06397 INTEGER :: mxcnt 06398 INTEGER :: nab 06399 INTEGER :: nkey 06400 INTEGER :: nkeys 06401 INTEGER :: nl 06402 INTEGER :: nmeth 06403 INTEGER :: npat 06404 INTEGER :: ntext 06405 INTEGER :: nums 06406 INTEGER :: matint 06407 06408 CHARACTER (LEN=*), INTENT(IN) :: text 06409 INTEGER, INTENT(IN) :: nline 06410 06411 PARAMETER (nkeys=9,nmeth=4) 06412 CHARACTER (LEN=16) :: methxt(nmeth) 06413 CHARACTER (LEN=16) :: keylst(nkeys) 06414 CHARACTER (LEN=32) :: keywrd 06415 CHARACTER (LEN=32) :: keystx 06416 DOUBLE PRECISION :: dnum(100) 06417 REAL :: plvs(3) ! vector array: real and ... 06418 INTEGER :: lpvs(3) ! ... integer 06419 EQUIVALENCE (plvs(1),lpvs(1)) 06420 06421 INTERFACE 06422 SUBROUTINE addItem(length,list,label,value) 06423 USE mpmod 06424 INTEGER, INTENT(IN OUT) :: length 06425 TYPE(listItem), DIMENSION(:), INTENT(IN OUT), ALLOCATABLE :: list 06426 INTEGER, INTENT(IN) :: label 06427 REAL, INTENT(IN) :: value 06428 END SUBROUTINE addItem 06429 END INTERFACE 06430 06431 DATA keylst/'unknown','parameter','constraint','wconstraint', & 06432 'measurement', 'method', & 06433 'mestimate', 'atleast','option'/ 06434 06435 !add number of iterations 06436 ! wconstraint like constraint 06437 ! measured r sigma label factor ...(more) 06438 06439 SAVE 06440 DATA methxt/'diagonalization','inversion','fullMINRES', 'sparseMINRES'/ 06441 DATA lkey/-1/ ! last keyword 06442 06443 ! ... 06444 nkey=-1 ! new keyword 06445 CALL rltext(text,ia,ib,nab) ! return indices for non-blank area 06446 IF(nab == 0) GOTO 10 06447 CALL ratext(text(1:nab),nums,dnum) ! translate text to DP numbers 06448 06449 IF(nums /= 0) nkey=0 06450 IF(keyb /= 0) THEN 06451 keywrd=text(keya:keyb) ! text is TEXT(KEYA)...TEXT(KEYB) 06452 ! WRITE(*,*) 'Keyword is ',KEYWRD 06453 06454 ! compare keywords 06455 06456 DO nkey=2,nkeys-1 ! loop over all pede keywords 06457 keystx=keylst(nkey) ! copy NKEY.th pede keyword 06458 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06459 IF(mat >= ntext-ntext/5) GO TO 10 06460 END DO 06461 06462 ! more comparisons 06463 06464 keystx='print' 06465 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06466 ! WRITE(*,*) KEYSTX,MAT,NTEXT 06467 ! IF(MAT.GE.NTEXT) THEN 06468 IF(mat >= (npat-npat/5)) THEN 06469 ! IF(MAT.GE.(NTEXT+NTEXT+3)/3) THEN 06470 mprint=1 06471 IF(nums > 0) mprint=IDNINT(dnum(1)) 06472 RETURN 06473 END IF 06474 06475 keystx='debug' 06476 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06477 IF(mat >= (npat-npat/5)) THEN 06478 ! IF(MAT.GE.(NTEXT+NTEXT+3)/3) THEN 06479 mdebug=3 06480 ! GF IF(NUMS.GT.0) MPRINT=DNUM(1) 06481 IF(nums > 0) mdebug=IDNINT(dnum(1)) 06482 IF(nums > 1) mdebg2=IDNINT(dnum(2)) 06483 RETURN 06484 END IF 06485 06486 keystx='entries' 06487 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06488 IF(mat >= (npat-npat/5)) THEN 06489 ! IF(MAT.GE.(NTEXT+NTEXT+3)/3) THEN 06490 mreqen=0 06491 IF(nums > 0) mreqen=IDNINT(dnum(1)) 06492 IF (mreqen <= 0) mreqen=1 06493 RETURN 06494 END IF 06495 06496 keystx='printrecord' 06497 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06498 IF(mat >= (npat-npat/5)) THEN 06499 IF(nums > 0) nrecpr=IDNINT(dnum(1)) 06500 IF(nums > 1) nrecp2=IDNINT(dnum(2)) 06501 RETURN 06502 END IF 06503 06504 keystx='maxrecord' 06505 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06506 IF(mat >= (npat-npat/5)) THEN 06507 IF (nums > 0.AND.dnum(1) > 0.) mxrec=IDNINT(dnum(1)) 06508 RETURN 06509 END IF 06510 06511 keystx='cache' 06512 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06513 IF(mat >= (npat-npat/5)) THEN 06514 IF (nums > 0.AND.dnum(1) >= 0.) ncache=IDNINT(dnum(1)) ! cache size, <0 keeps default 06515 IF (nums == 2.AND.dnum(2) > 0..AND.dnum(2) <= 1.0) & ! read cache fill level 06516 fcache(1)=SNGL(dnum(2)) 06517 IF (nums >= 4) THEN ! explicit cache splitting 06518 DO k=1,3 06519 fcache(k)=SNGL(dnum(k+1)) 06520 END DO 06521 END IF 06522 RETURN 06523 END IF 06524 06525 keystx='chisqcut' 06526 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06527 IF(mat >= (npat-npat/5)) THEN 06528 IF(nums == 0) THEN ! always 3-sigma cut 06529 chicut=1.0 06530 chirem=1.0 06531 ELSE 06532 chicut=SNGL(dnum(1)) 06533 IF(chicut < 1.0) chicut=-1.0 06534 IF(nums == 1) THEN 06535 chirem=1.0 ! 3-sigma cut, if not specified 06536 ELSE 06537 chirem=SNGL(dnum(2)) 06538 IF(chirem < 1.0) chirem=1.0 06539 IF(chicut >= 1.0) chirem=MIN(chirem,chicut) 06540 END IF 06541 END IF 06542 RETURN 06543 END IF 06544 06545 ! GF added: 06546 keystx='hugecut' 06547 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06548 IF(mat >= (npat-npat/5)) THEN 06549 IF(nums > 0) chhuge=SNGL(dnum(1)) 06550 IF(chhuge < 1.0) chhuge=1.0 ! at least (!!) 3-sigma 06551 RETURN 06552 END IF 06553 ! GF added end 06554 06555 keystx='localfit' 06556 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06557 IF(mat >= (npat-npat/5)) THEN 06558 IF(nums > 0) lfitnp=IDNINT(dnum(1)) 06559 IF(nums > 1) lfitbb=IDNINT(dnum(2)) 06560 RETURN 06561 END IF 06562 06563 keystx='regularization' 06564 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06565 IF(mat >= (npat-npat/5)) THEN 06566 nregul=1 06567 regula=SNGL(dnum(1)) 06568 IF(nums >= 2) regpre=SNGL(dnum(2)) 06569 RETURN 06570 END IF 06571 06572 keystx='regularisation' 06573 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06574 IF(mat >= (npat-npat/5)) THEN 06575 nregul=1 06576 regula=SNGL(dnum(1)) 06577 IF(nums >= 2) regpre=SNGL(dnum(2)) 06578 RETURN 06579 END IF 06580 06581 keystx='presigma' 06582 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06583 IF(mat >= (npat-npat/5)) THEN 06584 regpre=SNGL(dnum(1)) 06585 RETURN 06586 END IF 06587 06588 keystx='matiter' 06589 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06590 IF(mat >= (npat-npat/5)) THEN 06591 matrit=IDNINT(dnum(1)) 06592 RETURN 06593 END IF 06594 06595 keystx='matmoni' 06596 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06597 IF(mat >= (npat-npat/5)) THEN 06598 matmon=-1 06599 IF (nums > 0.AND.dnum(1) > 0.) matmon=IDNINT(dnum(1)) 06600 RETURN 06601 END IF 06602 06603 keystx='bandwidth' 06604 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06605 ! IF(MAT.GE.(NTEXT+NTEXT+3)/3) THEN 06606 IF(mat >= (npat-npat/5)) THEN 06607 IF(nums > 0) mbandw=IDNINT(dnum(1)) 06608 IF(mbandw < 0) mbandw=-1 06609 RETURN 06610 END IF 06611 06612 ! KEYSTX='outlierrejection' 06613 ! MAT=MATINT(TEXT(KEYA:KEYB),KEYSTX,NPAT,NTEXT) ! comparison 06614 ! WRITE(*,*) KEYSTX,MAT,(NTEXT+NTEXT)/3 06615 ! IF(MAT.GE.(NTEXT+NTEXT+NTEXT-2)/3) THEN 06616 ! IF(MAT.GE.(NPAT-NPAT/5)) THEN 06617 ! CHDFRJ=DNUM(1) 06618 ! IF(CHDFRJ.LT.3.0) CHDFRJ=100.0 06619 ! RETURN 06620 ! END IF 06621 06622 ! KEYSTX='outliersuppression' 06623 ! MAT=MATINT(TEXT(KEYA:KEYB),KEYSTX,NPAT,NTEXT) ! comparison 06624 ! WRITE(*,*) KEYSTX,MAT,(NTEXT+NTEXT)/3 06625 ! IF(MAT.GE.(NTEXT+NTEXT+NTEXT-2)/3) THEN 06626 ! IF(MAT.GE.(NPAT-NPAT/5)) THEN 06627 ! LHUBER=DNUM(1) 06628 ! IF(LHUBER.LE.2) LHUBER=2 ! at least 2 Huber iterations 06629 ! RETURN 06630 ! END IF 06631 06632 keystx='outlierdownweighting' 06633 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06634 ! WRITE(*,*) KEYSTX,MAT,(NTEXT+NTEXT)/3 06635 ! IF(MAT.GE.(NTEXT+NTEXT+NTEXT-2)/3) THEN 06636 IF(mat >= (npat-npat/5)) THEN 06637 lhuber=IDNINT(dnum(1)) 06638 IF(lhuber > 0.AND.lhuber <= 2) lhuber=2 ! at least 2 Huber iterations (if any) 06639 RETURN 06640 END IF 06641 06642 keystx='dwfractioncut' 06643 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06644 IF(mat >= (npat-npat/5)) THEN 06645 dwcut=SNGL(dnum(1)) 06646 IF(dwcut > 0.5) dwcut=0.5 06647 RETURN 06648 END IF 06649 06650 keystx='pullrange' 06651 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06652 IF(mat >= (npat-npat/5)) THEN 06653 prange=ABS(SNGL(dnum(1))) 06654 RETURN 06655 END IF 06656 06657 keystx='subito' 06658 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06659 IF(mat >= (npat-npat/5)) THEN 06660 isubit=1 06661 RETURN 06662 END IF 06663 06664 keystx='force' 06665 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06666 IF(mat >= (npat-npat/5)) THEN 06667 iforce=1 06668 RETURN 06669 END IF 06670 06671 keystx='memorydebug' 06672 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06673 IF(mat >= (npat-npat/5)) THEN 06674 memdbg=1 06675 IF (nums > 0.AND.dnum(1) > 0.0) memdbg=IDNINT(dnum(1)) 06676 RETURN 06677 END IF 06678 06679 keystx='globalcorr' 06680 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06681 IF(mat >= (npat-npat/5)) THEN 06682 igcorr=1 06683 RETURN 06684 END IF 06685 06686 keystx='threads' 06687 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06688 IF(mat >= (npat-npat/5)) THEN 06689 iomp=0 06690 !$ IOMP=1 06691 !$ IF (IOMP.GT.0) THEN 06692 !$ IF (NUMS.GE.1.AND.DNUM(1).GT.0.) MTHRD =IDNINT(DNUM(1)) 06693 !$ IF (NUMS.GE.2) THEN 06694 !$ MTHRDR=MTHRD 06695 !$ IF (DNUM(2).GT.0.) MTHRDR=IDNINT(DNUM(2)) 06696 !$ ENDIF 06697 !$ ELSE 06698 WRITE(*,*) 'WARNING: multithreading not available' 06699 !$ ENDIF 06700 RETURN 06701 END IF 06702 06703 keystx='compress' 06704 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06705 IF(mat >= (npat-npat/5)) THEN 06706 mcmprs=1 06707 RETURN 06708 END IF 06709 06710 keystx='errlabels' 06711 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06712 IF(mat >= (npat-npat/5).AND.mnrsel < 100) THEN 06713 nl=MIN(nums,100-mnrsel) 06714 DO k=1,nl 06715 lbmnrs(mnrsel+k)=IDNINT(dnum(k)) 06716 END DO 06717 mnrsel=mnrsel+nl 06718 RETURN 06719 END IF 06720 06721 keystx='pairentries' 06722 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06723 IF(mat >= (npat-npat/5)) THEN 06724 ! This option could be implemented to get rid of parameter pairs 06725 ! that have very few entries - to save matrix memory size. 06726 IF (nums > 0.AND.dnum(1) > 0.0) THEN 06727 mreqpe=IDNINT(dnum(1)) 06728 IF (nums >= 2.AND.dnum(2) >= dnum(1)) mhispe=IDNINT(dnum(2)) 06729 IF (nums >= 3.AND.dnum(3) > 0.0) msngpe=IDNINT(dnum(3)) 06730 icount=MAX(msngpe+1,mhispe) 06731 icount=MAX(mreqpe,icount) 06732 numbit=1 ! number of bits needed to count up to ICOUNT 06733 mxcnt=2 06734 DO k=1,30 06735 IF (icount >= mxcnt) THEN 06736 numbit=numbit+1 06737 mxcnt=mxcnt*2 06738 END IF 06739 END DO 06740 END IF 06741 RETURN 06742 END IF 06743 06744 keystx='wolfe' 06745 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06746 IF(mat >= (npat-npat/5)) THEN 06747 wolfc1=SNGL(dnum(1)) 06748 wolfc2=SNGL(dnum(2)) 06749 RETURN 06750 END IF 06751 06752 ! GF added: 06753 ! convergence tolerance for minres: 06754 keystx='mrestol' 06755 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06756 IF(mat >= (npat-npat/5)) THEN 06757 IF(nums > 0) THEN 06758 IF (dnum(1) < 1.0D-10.OR.dnum(1) > 1.0D-04) THEN 06759 WRITE(*,*) 'ERROR: need 1.0D-10 <= MRESTL ', & 06760 '<= 1.0D-04, but get ', dnum(1) 06761 ELSE 06762 mrestl=IDNINT(dnum(1)) 06763 END IF 06764 END IF 06765 RETURN 06766 END IF 06767 ! GF added end 06768 06769 keystx='nofeasiblestart' 06770 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06771 IF(mat >= (npat-npat/5)) THEN 06772 nofeas=1 ! do not make parameters feasible at start 06773 RETURN 06774 END IF 06775 06776 keystx='histprint' 06777 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06778 IF(mat >= ntext-ntext/10) THEN 06779 ! IF(MAT.GE.(NPAT-NPAT/5)) THEN 06780 nhistp=1 ! print histograms 06781 RETURN 06782 END IF 06783 06784 keystx='fortranfiles' 06785 mat=matint(text(ia:ib),keystx,npat,ntext) ! comparison 06786 IF(mat >= ntext-ntext/10) RETURN 06787 06788 keystx='Cfiles' 06789 mat=matint(text(ia:ib),keystx,npat,ntext) ! comparison 06790 IF(mat >= ntext-ntext/10) RETURN 06791 06792 keystx=keylst(1) 06793 nkey=1 ! unknown keyword 06794 IF(nums /= 0) nkey=0 06795 06796 WRITE(*,*) ' ' 06797 WRITE(*,*) '**************************************************' 06798 WRITE(*,*) ' ' 06799 WRITE(*,*) 'Unknown keyword(s): ',text(1:MIN(nab,50)) 06800 WRITE(*,*) ' ' 06801 WRITE(*,*) '**************************************************' 06802 WRITE(*,*) ' ' 06803 lunkno=lunkno+1 06804 06805 END IF 06806 ! result: NKEY = -1 blank 06807 ! NKEY = 0 numerical data, no text keyword or unknown 06808 ! NKEY > 0 keyword NKEY from list, keyword = KEYSTX 06809 06810 06811 ! content/lastcontent 06812 ! ------------------- 06813 ! blank -1 06814 ! data 0 06815 ! keyword 06816 ! unknown 1 06817 ! parameter 2 06818 ! constraint 3 06819 ! method 4 06820 ! atleast 5 06821 ! option 6 06822 ! cut 7 06823 06824 10 IF(nkey > 0) THEN ! new keyword 06825 lkey=nkey 06826 IF(lkey == 2) THEN ! parameter 06827 IF(nums == 3) THEN 06828 lpvs(1)=IDNINT(dnum(1)) ! label 06829 plvs(2)=SNGL(dnum(2)) ! start value 06830 plvs(3)=SNGL(dnum(3)) ! pre-sigma 06831 IF(lpvs(1) /= 0) THEN 06832 ! CALL megvec('7',plvs,3,0) ! always 3 words 06833 CALL addItem(lenParameters,listParameters,lpvs(1),plvs(2)) 06834 CALL addItem(lenPreSigmas,listPresigmas,lpvs(1),plvs(3)) 06835 ELSE 06836 WRITE(*,*) 'Line',nline,' error, label=',lpvs(1) 06837 END IF 06838 ELSE IF(nums /= 0) THEN 06839 kkey=1 ! switch to "unknown" ? 06840 WRITE(*,*) 'Wrong text in line',nline 06841 WRITE(*,*) 'Status: new parameter' 06842 WRITE(*,*) '> ',text(1:nab) 06843 END IF 06844 ELSE IF(lkey == 3.OR.lkey == 4) THEN ! constraint 06845 ! WRITE(*,*) 'Keyword is constraint!',NUMS,' numerical data' 06846 IF(nums >= 1.AND.nums <= 2) THEN ! start constraint 06847 ! WRITE(*,*) 'NUMs is ',NUMS 06848 lpvs(1)=0 06849 plvs(2)=SNGL(dnum(1)) ! r = r.h.s. value 06850 ! CALL megvec('8',plvs,2,0) 06851 CALL addItem(lenConstraints,listConstraints,lpvs(1),plvs(2)) 06852 ! WRITE(*,*) 'LPVS PLVS ',LPVS,PLVS 06853 lpvs(1)=-1 ! constraint 06854 IF(lkey == 4) lpvs(1)=-2 ! wconstraint (weighted) 06855 plvs(2)=0.0 06856 ! WRITE(*,*) 'LPVS PLVS ',LPVS,PLVS 06857 IF(nums == 2) plvs(2)=SNGL(dnum(2)) ! sigma 06858 ! CALL megvec('8',plvs,2,0) 06859 CALL addItem(lenConstraints,listConstraints,lpvs(1),plvs(2)) 06860 ! WRITE(*,*) 'LPVS PLVS ',LPVS,PLVS 06861 ELSE 06862 kkey=1 ! switch to "unknown" 06863 WRITE(*,*) 'Wrong text in line',nline 06864 WRITE(*,*) 'Status: new keyword (w)constraint' 06865 WRITE(*,*) '> ',text(1:nab) 06866 END IF 06867 ELSE IF(lkey == 5) THEN ! measurement 06868 IF(nums == 2) THEN ! start measurement 06869 lpvs(1)=0 06870 plvs(2)=SNGL(dnum(1)) ! r = r.h.s. value 06871 ! CALL megvec('9',plvs,2,0) 06872 CALL addItem(lenMeasurements,listMeasurements,lpvs(1),plvs(2)) 06873 lpvs(1)=-1 ! constraint 06874 plvs(2)=SNGL(dnum(2)) ! sigma 06875 ! CALL megvec('9',plvs,2,0) 06876 CALL addItem(lenMeasurements,listMeasurements,lpvs(1),plvs(2)) 06877 ELSE 06878 kkey=1 ! switch to "unknown" 06879 WRITE(*,*) 'Wrong text in line',nline 06880 WRITE(*,*) 'Status: new keyword measurement' 06881 WRITE(*,*) '> ',text(1:nab) 06882 END IF 06883 06884 ELSE IF(lkey == 6) THEN ! method 06885 miter=mitera 06886 IF(nums >= 1) miter=IDNINT(dnum(1)) 06887 IF(miter >= 1) mitera=miter 06888 dflim=SNGL(dnum(2)) 06889 lkey=0 06890 DO i=1,nmeth 06891 keystx=methxt(i) 06892 mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison 06893 IF(mat >= ntext-ntext/5) THEN 06894 IF(i == 1) THEN ! diagonalization 06895 metsol=2 06896 matsto=1 06897 ELSE IF(i == 2) THEN ! inversion 06898 metsol=1 06899 matsto=1 06900 ELSE IF(i == 3) THEN ! fullMINRES 06901 metsol=3 06902 matsto=1 06903 ! METSOL=4 !!! 06904 ELSE IF(i == 4) THEN ! sparseMINRES 06905 ! METSOL=4 06906 metsol=3 06907 matsto=2 06908 06909 END IF 06910 END IF 06911 END DO 06912 END IF 06913 ELSE IF(nkey == 0) THEN ! data for continuation 06914 IF(lkey == 2) THEN ! parameter 06915 IF(nums >= 3) THEN ! store data from this line 06916 lpvs(1)=IDNINT(dnum(1)) ! label 06917 plvs(2)=SNGL(dnum(2)) ! start value 06918 plvs(3)=SNGL(dnum(3)) ! pre-sigma 06919 IF(lpvs(1) /= 0) THEN 06920 ! CALL megvec('7',plvs,3,0) ! always 3 words 06921 CALL addItem(lenParameters,listParameters,lpvs(1),plvs(2)) 06922 CALL addItem(lenPreSigmas,listPresigmas,lpvs(1),plvs(3)) 06923 ELSE 06924 WRITE(*,*) 'Line',nline,' error, label=',lpvs(1) 06925 END IF 06926 ELSE IF(nums > 1.AND.nums < 3) THEN 06927 kkey=1 ! switch to "unknown" ? 06928 WRITE(*,*) 'Wrong text in line',nline 06929 WRITE(*,*) 'Status continuation parameter' 06930 WRITE(*,*) '> ',text(1:nab) 06931 END IF 06932 06933 ELSE IF(lkey == 3.OR.lkey == 4) THEN ! (w)constraint 06934 ier=0 06935 DO i=1,nums,2 06936 label=IDNINT(dnum(i)) 06937 IF(label <= 0) ier=1 06938 END DO 06939 IF(MOD(nums,2) /= 0) ier=1 ! reject odd number 06940 IF(ier == 0) THEN 06941 DO i=1,nums,2 06942 lpvs(1)=IDNINT(dnum(i)) ! label 06943 plvs(2)=SNGL(dnum(i+1)) ! factor 06944 ! CALL megvec('8',plvs,2,0) 06945 CALL addItem(lenConstraints,listConstraints,lpvs(1),plvs(2)) 06946 END DO 06947 ELSE 06948 kkey=0 06949 WRITE(*,*) 'Wrong text in line',nline 06950 WRITE(*,*) 'Status continuation (w)constraint' 06951 WRITE(*,*) '> ',text(1:nab) 06952 END IF 06953 06954 ELSE IF(lkey == 5) THEN ! measurement 06955 ! WRITE(*,*) 'continuation < ',NUMS 06956 ier=0 06957 DO i=1,nums,2 06958 label=IDNINT(dnum(i)) 06959 IF(label <= 0) ier=1 06960 END DO 06961 IF(MOD(nums,2) /= 0) ier=1 ! reject odd number 06962 ! WRITE(*,*) 'IER NUMS ',IER,NUMS 06963 IF(ier == 0) THEN 06964 DO i=1,nums,2 06965 lpvs(1)=IDNINT(dnum(i)) ! label 06966 plvs(2)=SNGL(dnum(i+1)) ! factor 06967 ! CALL megvec('9',plvs,2,0) 06968 CALL addItem(lenMeasurements,listMeasurements,lpvs(1),plvs(2)) 06969 END DO 06970 ELSE 06971 kkey=0 06972 WRITE(*,*) 'Wrong text in line',nline 06973 WRITE(*,*) 'Status continuation measurement' 06974 WRITE(*,*) '> ',text(1:nab) 06975 END IF 06976 06977 END IF 06978 END IF 06979 END SUBROUTINE intext 06980 06988 SUBROUTINE addItem(length,list,label,value) 06989 USE mpdalc 06990 06991 INTEGER, INTENT(IN OUT) :: length 06992 TYPE(listItem), DIMENSION(:), INTENT(IN OUT), ALLOCATABLE :: list 06993 INTEGER, INTENT(IN) :: label 06994 REAL, INTENT(IN) :: value 06995 06996 INTEGER(kind=large) :: newSize 06997 INTEGER(kind=large) :: oldSize 06998 TYPE(listItem), DIMENSION(:), ALLOCATABLE :: tempList 06999 07000 IF (label > 0.AND.value == 0.) RETURN ! skip zero for valid labels 07001 IF (length == 0 ) THEN ! initial list with size = 100 07002 newSize = 100 07003 CALL mpalloc(list,newSize,' list ') 07004 ENDIF 07005 oldSize=size(list,kind=large) 07006 IF (length >= oldSize) THEN ! increase sizeby 20% + 100 07007 newSize = oldSize + oldSize/5 + 100 07008 CALL mpalloc(tempList,oldSize,' temp. list ') 07009 tempList=list 07010 CALL mpdealloc(list) 07011 CALL mpalloc(list,newSize,' list ') 07012 list(1:oldSize)=tempList(1:oldSize) 07013 CALL mpdealloc(tempList) 07014 ENDIF 07015 ! add to end of list 07016 length=length+1 07017 list(length)%label=label 07018 list(length)%value=value 07019 07020 END SUBROUTINE addItem 07021 07023 SUBROUTINE mstart(text) 07024 USE mpmod, ONLY: textl 07025 07026 IMPLICIT NONE 07027 INTEGER :: i 07028 INTEGER :: ka 07029 INTEGER :: kb 07030 INTEGER :: l 07031 CHARACTER (LEN=*), INTENT(IN) :: text 07032 CHARACTER (LEN=16) :: textc 07033 SAVE 07034 ! ... 07035 DO i=1,74 07036 textl(i:i)='_' 07037 END DO 07038 l=LEN(text) 07039 ka=(74-l)/2 07040 kb=ka+l-1 07041 textl(ka:kb)=text(1:l) 07042 WRITE(*,*) ' ' 07043 WRITE(*,*) textl 07044 WRITE(*,*) ' ' 07045 textc=text(1:l)//'-end' 07046 07047 DO i=1,74 07048 textl(i:i)='_' 07049 END DO 07050 l=l+4 07051 ka=(74-l)/2 07052 kb=ka+l-1 07053 textl(ka:kb)=textc(1:l) 07054 RETURN 07055 END SUBROUTINE mstart 07056 07058 SUBROUTINE mend 07059 USE mpmod, ONLY: textl 07060 07061 IMPLICIT NONE 07062 WRITE(*,*) ' ' 07063 WRITE(*,*) textl 07064 CALL petime 07065 WRITE(*,*) ' ' 07066 END SUBROUTINE mend 07067 07074 07075 SUBROUTINE mvopen(lun,fname) 07076 IMPLICIT NONE 07077 INTEGER :: l 07078 INTEGER, INTENT(IN) :: lun 07079 CHARACTER (LEN=*), INTENT(IN) :: fname 07080 CHARACTER (LEN=33) :: nafile 07081 CHARACTER (LEN=33) :: nbfile 07082 LOGICAL :: ex 07083 SAVE 07084 ! ... 07085 l=LEN(fname) 07086 IF(l > 32) STOP 'File name too long ' 07087 nafile=fname 07088 nafile(l+1:l+1)='~' 07089 07090 INQUIRE(FILE=nafile(1:l),EXIST=ex) 07091 IF(ex) THEN 07092 INQUIRE(FILE=nafile(1:l+1),EXIST=ex) 07093 IF(ex) THEN 07094 CALL system('rm '//nafile) 07095 END IF 07096 nbfile=nafile 07097 nafile(l+1:l+1)=' ' 07098 CALL system('mv '//nafile//nbfile) 07099 END IF 07100 OPEN(UNIT=lun,FILE=fname) 07101 END SUBROUTINE mvopen 07102 07106 07107 SUBROUTINE petime 07108 IMPLICIT NONE 07109 REAL :: ta(2) 07110 REAL :: rst 07111 REAL :: delta 07112 REAL :: rstp 07113 REAL :: secnd1 07114 REAL :: secnd2 07115 INTEGER :: ncount 07116 INTEGER :: nhour1 07117 INTEGER :: minut1 07118 INTEGER :: nsecd1 07119 INTEGER :: nhour2 07120 INTEGER :: minut2 07121 INTEGER :: nsecd2 07122 07123 SAVE 07124 DATA ncount/0/ 07125 ! ... 07126 ncount=ncount+1 07127 CALL etime(ta,rst) 07128 IF(ncount > 1) THEN 07129 delta=rst 07130 nsecd1=INT(delta) ! -> integer 07131 nhour1=nsecd1/3600 07132 minut1=nsecd1/60-60*nhour1 07133 secnd1=delta-60*(minut1+60*nhour1) 07134 delta=rst-rstp 07135 nsecd2=INT(delta) ! -> integer 07136 nhour2=nsecd2/3600 07137 minut2=nsecd2/60-60*nhour2 07138 secnd2=delta-60*(minut2+60*nhour2) 07139 WRITE(*,101) nhour1,minut1,secnd1, nhour2,minut2,secnd2 07140 END IF 07141 07142 rstp=rst 07143 RETURN 07144 101 FORMAT(i4,' h',i3,' min',f5.1,' sec total',18X,'elapsed', & 07145 i4,' h',i3,' min',f5.1,' sec') 07146 END SUBROUTINE petime ! print 07147 07148 ! ----- accurate summation ----(from mpnum) --------------------------------- 07149 07153 07154 SUBROUTINE addsum(add) 07155 USE mpmod 07156 07157 IMPLICIT NONE 07158 DOUBLE PRECISION:: add 07159 INTEGER::nadd 07160 ! ... 07161 nadd=INT(add) ! convert to integer 07162 accurateNsum=accurateNsum+nadd ! sum integer 07163 accurateDsum=accurateDsum+(add-dfloat(nadd)) ! sum remainder 07164 IF(accurateDsum > 16.0D0) THEN ! + - 16 07165 accurateDsum=accurateDsum-16.0D0 07166 accurateNsum=accurateNsum+16 07167 END IF 07168 IF(accurateNsum > nexp20) THEN ! if > 2^20: + - 2^20 07169 accurateNexp=accurateNexp+1 07170 accurateNsum=accurateNsum-nexp20 07171 END IF 07172 RETURN 07173 END SUBROUTINE addsum 07174 07178 07179 SUBROUTINE getsum(asum) 07180 USE mpmod 07181 07182 IMPLICIT NONE 07183 DOUBLE PRECISION, INTENT(OUT) ::asum 07184 asum=(accurateDsum+dfloat(accurateNsum))+dfloat(accurateNexp)*dfloat(nexp20) 07185 accurateDsum=0.0D0 07186 accurateNsum=0 07187 accurateNexp=0 07188 RETURN 07189 END SUBROUTINE getsum 07190
1.8.0
