C++********************************************************************* C C DSGR_P.F C CROSRNG_E SPEEDS UP AUG 04 ARDEAN LEITH C AP_END CALL HAS PARLIST OCT 04 ARDEAN LEITH C SPIDER REF_RINGS FILE FEB 05 ARDEAN LEITH C CCROT STATISTICS FEB 05 ARDEAN LEITH C NON-INCORE PARALLEL FEB 05 ARDEAN LEITH C APRINGS_ONE MAR 05 ARDEAN LEITH C FFTW_PLANS MAR 08 ARDEAN LEITH C APRINGS_ONE_NEW MAR 08 ARDEAN LEITH C MPI PARTAB BUG FIXED OCT 08 ARDEAN LEITH C C ********************************************************************** C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM. AUTHOR: J.FRANK * C=* Copyright (C) 1985-2008 Health Research Inc. * C=* * C=* HEALTH RESEARCH INCORPORATED (HRI), * C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455. * C=* * C=* Email: spider@wadsworth.org * C=* * C=* This program is free software; you can redistribute it and/or * C=* modify it under the terms of the GNU General Public License as * C=* published by the Free Software Foundation; either version 2 of the * C=* License, or (at your option) any later version. * C=* * C=* This program is distributed in the hope that it will be useful, * C=* but WITHOUT ANY WARRANTY; without even the implied warranty of * C=* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * C=* General Public License for more details. * C=* * C=* You should have received a copy of the GNU General Public License * C=* along with this program; if not, write to the * C=* Free Software Foundation, Inc., * C=* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * C=* * C ********************************************************************** C C DSGR_P(ILIST,NIMA,ILIP,NIDI, C NSAM,NROW,NR,LENTT,RANGE, C NRING,LCIRC,NUMR,CIRCREF,CIRCREF_IN_CORE, C MODE,REFANGDOC,EXPANGDOC,SCRFILE,FFTW_PLANS, C REFPAT,EXPPAT,CKMIRROR,CTYPE,ISHRANGE,LUNDOC) C C PURPOSE: FIND ROTATIONAL AND SHIFT PARAMETERS TO ALIGN A SERIES OF C REFERENCE IMAGES WITH SAMPLE IMAGES C C ORDER OF PROCESSING: C 1 - CONVERT EACH REFERENCE IMAGE TO RINGS, DO THE FFTS C FOR ALL THE RINGS, APPLY WEIGHTS TO THE RINGS, STORE IT IN CIRCREF. C IN ADDITION, HIGHEST FREQUENCY FOR ALL THE RINGS EXCEPT C MAXRING ARE DIVIDED BY 2. C 2 - CONVERT EACH INPUT IMAGE TO RINGS, DO THE FFTS, C COMPARE EACH INPUT IMAGE WITH ALL THE REFERENCE IMAGES C USING CROSRNG_**. C SINCE Y WERE PRE-WEIGHTED THE RESULT ARE ALREADY CORRECT. C C PARAMETERS: C ILIST LIST OF REF. IMAGE FILE NUMBERS (INPUT) C NIMA NO. OF IMAGES (INPUT) C ILIP LIST OF EXP. IMAGE FILE NUMBERS (INPUT) C NIDI NO. OF IMAGES (INPUT) C REFANGDOC REF. ANGLES FILE NAME (INPUT) C EXPANGDOC EXP. ANGLES FILE NAME (INPUT) C REFPAT REF. IMAGE SERIES FILE TEMPLATE (INPUT) C EXPPAT EXP. IMAGE SERIES FILE TEMPLATE (INPUT) C C OPERATIONS: 'AP REF', 'AP RD', 'AP RN' C C--********************************************************************* #ifndef USE_MPI C --------------------- NON-MPI CODE -------------------------------- SUBROUTINE DSGR_P(ILIST,NIMA, ILIP,NIDI, & NSAM,NROW,LENTT,RANGE,ANGDIFTHR, & NRING,LCIRC,NUMR,CIRCREF,CIRCREF_IN_CORE, & MODE,REFANGDOC,EXPANGDOC,SCRFILE,FFTW_PLANS, & REFPAT,EXPPAT,CKMIRROR,CTYPE,ISHRANGE,LUNDOC) INCLUDE 'CMLIMIT.INC' INCLUDE 'CMBLOCK.INC' INTEGER, DIMENSION(NIMA) :: ILIST INTEGER, DIMENSION(NIDI) :: ILIP INTEGER, DIMENSION(3,NRING) :: NUMR REAL, DIMENSION(LCIRC,NIMA) :: CIRCREF CHARACTER(LEN=1) :: MODE CHARACTER (LEN=*) :: REFANGDOC CHARACTER (LEN=*) :: EXPANGDOC CHARACTER (LEN=*) :: SCRFILE CHARACTER (LEN=*) :: REFPAT,EXPPAT CHARACTER (LEN=*) :: CTYPE CHARACTER (LEN=74) :: COMMENT DOUBLE PRECISION :: CCROTD CHARACTER(LEN=1) :: NULL LOGICAL :: CIRCREF_IN_CORE LOGICAL :: CKMIRROR LOGICAL :: MIRRORNEW LOGICAL :: GOTREFANG LOGICAL :: ONLYMIRROR LOGICAL :: LIMITRANGE LOGICAL :: MIRRORED INTEGER *8 :: FFTW_PLANS(*) C ALLOCATABLE ARRAYS REAL, ALLOCATABLE, DIMENSION(:) :: XBUF C AUTOMATIC ARRAYS DOUBLE PRECISION, DIMENSION(NIMA) :: TOTMIN,TOTMIR DOUBLE PRECISION, DIMENSION(LENTT) :: TT REAL, DIMENSION(NIMA) :: ROTANGT,ROTANGM REAL, DIMENSION(LCIRC) :: CIRCEXP REAL, DIMENSION(6) :: DLIST REAL, DIMENSION(3) :: EXPDIR REAL, DIMENSION(8,NIDI) :: ANGEXP REAL, DIMENSION(3,NIMA) :: ANGREF,REFDIR INTEGER, DIMENSION(NIMA) :: LCG INTEGER :: NSAID = 0 INTEGER, PARAMETER :: NLISTMAX = 15 REAL, DIMENSION(NLISTMAX) :: PARLIST PARAMETER (QUADPI = 3.1415926535897932384626) PARAMETER (DGR_TO_RAD = (QUADPI/180)) DATA INPIC,LUNANG,LUNRING/77,78,50/ MYPID = -1 ! NOT USING MPI NULL = CHAR(0) MAXRIN = NUMR(3,NRING) ! ACTUAL LENGTH OF LONGEST RING #ifdef SP_LIBFFTW3 MAXRIN = NUMR(3,NRING) - 2 #endif C THIS ALTERS RANGE! RANGE = COS(RANGE*DGR_TO_RAD) C FIND NUMBER OF OMP THREADS CALL GETTHREADS(NUMTH) IF (MODE .EQ. 'H') THEN DIVAS = 180.0 ELSE DIVAS = 360.0 ENDIF C MAKE XBUF AT LEAST NSAM*NROW FOR USE BY AP_SHIFT MWANTX = NSAM * NROW ALLOCATE(XBUF(MWANTX),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'XBUF',MWANTX) GOTO 9999 ENDIF C LOAD REF. PROJ. ANGLES (ANGREF) FROM DOC. FILE (REFANGDOC) OR C REF. IMAGE FILE (REFPAT) HEAD CALL AP_GETANGA(ILIST,NIMA,0,REFANGDOC,REFPAT, & INPIC,LUNANG,3,ANGREF,NGOTREF,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C CONVERT REF. ANGLES TO UNITARY DIRECTIONAL VECTORS (REFDIR). CALL AP_GETSATA(ANGREF,REFDIR,3,NIMA,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 GOTREFANG = .TRUE. C READ REFERENCE IMAGES INTO REFERENCE RINGS (CIRCREF) ARRAY OR C CREATE REFERENCE RINGS FILE FOR LATER READING CALL APRINGS_NEW(ILIST,NIMA, NSAM,NROW, & NRING,LCIRC,NUMR,MODE,FFTW_PLANS, & REFPAT,INPIC,CIRCREF,CIRCREF_IN_CORE, & LUNRING,SCRFILE,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 IF (NSAID .LE. 0) THEN IF (CIRCREF_IN_CORE) THEN WRITE(NOUT,91)NUMTH 91 FORMAT(' Ref. rings in core, Threads: ',I4) ELSE WRITE(NOUT,92)NUMTH 92 FORMAT(' Ref. rings not in core, Threads: ',I4) ENDIF NSAID = NSAID + 1 ENDIF C LOAD EXP. PROJ. ANGLES & ALIGNMENT PARAMETERS (ANGEXP) C FROM DOC. FILE (EXPANGDOC) NGOTPAR = 0 IF (EXPANGDOC .NE. NULL) THEN C LOAD EXP. PROJ. ANGLES & ALIGNMENT PARAMETERS (ANGEXP) C FROM DOC. FILE (EXPANGDOC) CALL AP_GETANGA(ILIP,NIDI,0,EXPANGDOC,EXPPAT, & INPIC,LUNANG,8,ANGEXP,NGOTPAR,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ELSE ANGEXP = 0.0 ENDIF C INITIALIZE CCROT STATISTICS CCROTAVG = 0.0 IMPROVCCROT = 0 CCROTIMPROV = 0.0 IWORSECCROT = 0 CCROTWORSE = 0.0 ANGDIFAVG = 0.0 IBIGANGDIF = 0 C CALCULATE DIMENSIONS FOR APRINGS CNS2 = NSAM / 2 + 1 CNR2 = NROW / 2 + 1 DO IEXP=1,NIDI C LOOP OVER ALL EXPERIMENTAL (SAMPLE) IMAGES C CONVERT EXP. ANGLE TO UNITARY DIRECTIONAL VECTORS (EXPDIR). CALL AP_GETSATA(ANGEXP(1,IEXP),EXPDIR,8,1,IRTFLG) C LOAD CURRENT EXPERIMENTAL IMAGE INTO ARRAY XBUF CALL AP_GETDATS(ILIP,NIDI,NSAM,NROW, & 1,EXPPAT,INPIC, IEXP,IEXP, & XBUF, IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C EXTRACT EXP. IMAGE POLAR COORD. RINGS, NORMALIZE & FFT THEM CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, XBUF,.FALSE., & MODE,NUMR,NRING,LCIRC, 0.0,FFTW_PLANS, & CIRCEXP,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C DETERMINE WHICH REF IMAGES ARE TO BE COMPARED NIMALCG = 0 IEND = NIMA CCROTD = -1.0D20 LIMITRANGE = RANGE .LT. 1.0 IF (LIMITRANGE) THEN C DETERMINE WHICH REF IMAGES ARE TO BE COMPARED DO IMI=1,NIMA C LOOP OVER ALL REF. IMAGES C ABS - DIRECTIONS AT 180 DEGREES ARE DIFFERENT C (- DO NOT CHECK MIRRORED) C DT NEAR 1.0 = NOT-MIRRORED, DT NEAR -1.0 = MIRRORED DT = (EXPDIR(1) * REFDIR(1,IMI) + & EXPDIR(2) * REFDIR(2,IMI) + & EXPDIR(3) * REFDIR(3,IMI)) DTABS = ABS(DT) IF (DTABS .GE. RANGE) THEN C NON-MIRRORED OR MIRRORED REF. IS WITHIN RANGE NIMALCG = NIMALCG + 1 LCG(NIMALCG) = IMI IF (DT .LT. 0) LCG(NIMALCG) = -IMI ENDIF ENDDO IF (NIMALCG .LE. 0) THEN C NO REF. IMAGE WITHIN COMPARISON ANGLE IREF = 0 GOTO 1000 ENDIF IEND = NIMALCG ENDIF C REF. IMAGES FOUND WITHIN COMPARISION RANGE IF (CIRCREF_IN_CORE) THEN C USE CIRCREF FOR REFERENCE RINGS c$omp parallel do private(imil,imi,mirrored) DO IMIL=1,IEND IMI = IMIL IF (LIMITRANGE) IMI = ABS(LCG(IMIL)) IF ((CKMIRROR .AND. LIMITRANGE) .OR. & (.NOT. CKMIRROR)) THEN IF (.NOT. CKMIRROR) MIRRORED = .FALSE. IF (LIMITRANGE) MIRRORED = (LCG(IMIL) .LT. 0) C CHECK EITHER MIRRORED OR NON-MIRRORED POSITION CALL CROSRNG_EP_NEW(CIRCREF(1,IMI),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TT,MIRRORED,FFTW_PLANS(1)) ELSE C CHECK BOTH NON-MIRRORED & MIRRORED POSITIONS CALL CROSRNG_MSP_NEW(CIRCREF(1,IMI),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TOTMIR(IMI),ROTANGM(IMI), & TT,FFTW_PLANS(1)) ENDIF ENDDO c$omp end parallel do ELSE C USE REFERENCE RINGS FILE (MIGHT BE AN INCORE FILE) c$omp parallel do private(imil,imi,mirrored,ithread) c$omp& schedule(static,1) DO IMIL=1,IEND !LOOP OVER ALL REFERENCE IMAGES IMI = IMIL IF (LIMITRANGE) IMI = ABS(LCG(IMIL)) C FIND THREAD NUMBER ITHREAD = MOD((IMIL-1),NUMTH) + 1 C FILL CIRCREF FROM REFERENCE RINGS FILE c$omp critical CALL REDLIN(LUNRING,CIRCREF(1,ITHREAD),LCIRC,IMI) c$omp end critical IF ((CKMIRROR .AND. LIMITRANGE) .OR. & (.NOT. CKMIRROR)) THEN IF (.NOT. CKMIRROR) MIRRORED = .FALSE. IF (LIMITRANGE) MIRRORED = (LCG(IMIL) .LT. 0) C CHECK EITHER MIRRORED OR NON-MIRRORED POSITION CALL CROSRNG_EP_NEW(CIRCREF(1,ITHREAD),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TT,MIRRORED,FFTW_PLANS(1)) ELSE C CHECK BOTH NON-MIRRORED & MIRRORED POSITIONS CALL CROSRNG_MSP_NEW(CIRCREF(1,ITHREAD),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TOTMIR(IMI),ROTANGM(IMI), & TT,FFTW_PLANS(1)) ENDIF ENDDO comp end parallel do ENDIF C LOOP OVER ALL RELEVANT REF. IMAGES IREF = 0 CCROTD = -1.0D20 DO IMIL=1,IEND IMI = IMIL IF (NIMALCG .GT. 0) IMI = ABS(LCG(IMIL)) IF (TOTMIN(IMI) .GE. CCROTD) THEN C GOOD MATCH WITH TOTMIN (MIRRORED OR NOT) POSITION CCROTD = TOTMIN(IMI) RANGNEW = ROTANGT(IMI) MIRRORNEW = (LIMITRANGE .AND. (LCG(IMIL) .LT. 0)) IREF = IMI ENDIF IF (CKMIRROR .AND. .NOT. LIMITRANGE) THEN C HAVE TO COMPARE WITH MIRRORED POSITION IF (TOTMIR(IMI) .GE. CCROTD) THEN C GOOD MATCH, MIRRORED POSITION IS BETTER CCROTD = TOTMIR(IMI) RANGNEW = ROTANGM(IMI) MIRRORNEW = .TRUE. IREF = IMI ENDIF ENDIF ! END OF: IF (CKMIRROR) ENDDO ! END OF: DO IMIL=1,IEND 1000 RANGNEW = (RANGNEW-1) / MAXRIN * DIVAS CCROT = CCROTD IMGEXP = ILIP(IEXP) PEAKV = 0.0 XSHNEW = 0.0 YSHNEW = 0.0 IF (IREF .LE. 0) THEN C NO NEARBY REFERENCE IMAGE IMGREF = 0 C IREFT IS FOR REFDIR INDEX IREFT = 1 ELSE IMGREF = ILIST(IREF) C IREFT IS FOR REFDIR INDEX IREFT = IREF ENDIF C AP_END DETERMINES SHIFT PARAMETERS AND WRITES TO DOC FILE NPROJ = NIMA IF (LIMITRANGE) NPROJ = NIMALCG CALL AP_END(IEXP,IMGEXP,IMGREF, & ANGREF(1,IREFT),REFDIR(1,IREFT), & ANGEXP(1,IEXP),EXPDIR,ISHRANGE, & GOTREFANG, NGOTPAR,NSAM,NROW,CCROT,PEAKV, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW,EXPPAT,REFPAT, & NPROJ, CTYPE, XBUF,LUNDOC,PARLIST) CCROTAVG = CCROTAVG + CCROT IF (NGOTPAR .GE. 8) THEN C COMPILE CCROT CHANGE STATISTICS ANGDIF = PARLIST(10) IF (ANGDIF .GT. ANGDIFTHR)IBIGANGDIF = IBIGANGDIF + 1 CCROTLAS = ANGEXP(8,IEXP) ANGDIFAVG = ANGDIFAVG + PARLIST(10) CCROTAVG = CCROTAVG + CCROT IF (CCROT .GE. CCROTLAS) THEN IMPROVCCROT = IMPROVCCROT + 1 CCROTIMPROV = CCROTIMPROV + CCROT ELSE IWORSECCROT = IWORSECCROT + 1 CCROTWORSE = CCROTWORSE + CCROT ENDIF ENDIF ENDDO IF (NIDI .GT. 1) THEN C SAVE CCROT & ANGULAR DISPLACEMENT STATISTICS CALL AP_STAT(NIDI,ANGDIFTHR,IBIGANGDIF, & ANGDIFAVG, CCROTAVG, & IMPROVCCROT,CCROTIMPROV, & IWORSECCROT,CCROTWORSE,LUNDOC) ENDIF 9999 CLOSE(LUNANG) IF (.NOT. CIRCREF_IN_CORE) CLOSE(LUNRING) IF (ALLOCATED(XBUF)) DEALLOCATE(XBUF) END #else C------------------------ MPI SPECIFIC SUBROUTINE -------------------- SUBROUTINE DSGR_P(ILIST,NIMA,ILIP,NIDI, & NSAM,NROW,LENTT,RANGE,ANGDIFTHR, & NRING,LCIRC,NUMR,CIRCREF,CIRCREF_IN_CORE, & MODE,REFANGDOC,EXPANGDOC,SCRFILE,FFTW_PLANS, & REFPAT,EXPPAT,CKMIRROR,CTYPE,ISHRANGE,LUNDOC) INCLUDE 'CMLIMIT.INC' INCLUDE 'CMBLOCK.INC' INTEGER, DIMENSION(NIMA) :: ILIST INTEGER, DIMENSION(NIDI) :: ILIP INTEGER, DIMENSION(3,NRING) :: NUMR REAL, DIMENSION(LCIRC,NIMA) :: CIRCREF CHARACTER(LEN=1) :: MODE CHARACTER (LEN=*) :: REFANGDOC CHARACTER (LEN=*) :: EXPANGDOC CHARACTER (LEN=*) :: SCRFILE CHARACTER (LEN=*) :: REFPAT,EXPPAT CHARACTER (LEN=*) :: CTYPE CHARACTER (LEN=74) :: COMMENT DOUBLE PRECISION :: CCROTD CHARACTER(LEN=1) :: NULL LOGICAL :: CIRCREF_IN_CORE LOGICAL :: CKMIRROR LOGICAL :: MIRRORNEW LOGICAL :: GOTREFANG LOGICAL :: ONLYMIRROR LOGICAL :: LIMITRANGE LOGICAL :: MIRRORED INTEGER *8 :: FFTW_PLANS(*) C ALLOCATABLE ARRAYS REAL, ALLOCATABLE, DIMENSION(:) :: XBUF C AUTOMATIC ARRAYS DOUBLE PRECISION, DIMENSION(NIMA) :: TOTMIN,TOTMIR DOUBLE PRECISION, DIMENSION(LENTT) :: TT REAL, DIMENSION(NIMA) :: ROTANGT,ROTANGM REAL, DIMENSION(LCIRC) :: CIRCEXP REAL, DIMENSION(6) :: DLIST REAL, DIMENSION(3) :: EXPDIR REAL, DIMENSION(8,NIDI) :: ANGEXP REAL, DIMENSION(3,NIMA) :: ANGREF,REFDIR INTEGER, DIMENSION(NIMA) :: LCG INTEGER :: NSAID = 0 INTEGER, PARAMETER :: NLISTMAX = 15 REAL, DIMENSION(NLISTMAX) :: PARLIST PARAMETER (QUADPI = 3.1415926535897932384626) PARAMETER (DGR_TO_RAD = (QUADPI/180)) DATA INPIC,LUNANG,LUNRING/77,78,50/ INCLUDE 'mpif.h' INTEGER :: IGRPS(3) INTEGER, ALLOCATABLE, DIMENSION(:) :: PSIZE INTEGER, ALLOCATABLE, DIMENSION(:) :: NBASE REAL, ALLOCATABLE, DIMENSION(:,:) :: PARTAB,PARTABLOC REAL, ALLOCATABLE, DIMENSION(:,:) :: ALOC, ABUF #ifdef MPI_DEBUG DOUBLE PRECISION TCOM0, TCOM1 #endif LOGICAL :: USEBCAST COMMON /COMM_MPI/USEBCAST ICOMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(ICOMM, MYPID, MPIERR) CALL MPI_COMM_SIZE(ICOMM, NPROCS, MPIERR) NULL = CHAR(0) NR = NUMR(1,NRING) ! OUTER RING NUMBER MAXRIN = NUMR(3,NRING) ! ACTUAL LENGTH OF LONGEST RING #ifdef SP_LIBFFTW3 MAXRIN = NUMR(3,NRING) - 2 #endif C THIS ALTERS RANGE! RANGE = COS(RANGE*DGR_TO_RAD) C FIND NUMBER OF OMP THREADS CALL GETTHREADS(NUMTH) IF (MODE .EQ. 'H') THEN DIVAS = 180.0 ELSE DIVAS = 360.0 ENDIF C MAKE XBUF AT LEAST NSAM*NROW FOR USE BY AP_SHIFTS MWANTX = NSAM*NROW ALLOCATE(XBUF(MWANTX),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'XBUF',MWANTX) GOTO 9999 ENDIF C LOAD REF. PROJ. ANGLES (ANGREF) FROM DOC. FILE (REFANGDOC) OR C REF. IMAGE FILE (REFPAT) HEAD CALL AP_GETANGA(ILIST,NIMA,0,REFANGDOC,REFPAT, & INPIC,LUNANG,3,ANGREF,NGOTREF,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C CONVERT REF. ANGLES TO UNITARY DIRECTIONAL VECTORS (REFDIR). CALL AP_GETSATA(ANGREF,REFDIR,3,NIMA,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 GOTREFANG = .TRUE. C READ REFERENCE IMAGES INTO REFERENCE RINGS (CIRCREF) ARRAY OR C CREATE REFERENCE RINGS FILE FOR LATER READING CALL APRINGS_NEW(ILIST,NIMA, NSAM,NROW, & NRING,LCIRC,NUMR,MODE,FFTW_PLANS, & REFPAT,INPIC,CIRCREF,CIRCREF_IN_CORE, & LUNRING,SCRFILE,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 IF (NSAID .LE. 0) THEN IF (CIRCREF_IN_CORE) THEN WRITE(NOUT,91)NUMTH 91 FORMAT(' Ref. rings in core, Threads: ',I4) ELSE WRITE(NOUT,92)NUMTH 92 FORMAT(' Ref. rings not in core, Threads: ',I4) ENDIF NSAID = NSAID + 1 ENDIF C LOAD EXP. PROJ. ANGLES & ALIGNMENT PARAMETERS (ANGEXP) C FROM DOC. FILE (EXPANGDOC) NGOTPAR = 0 IF (EXPANGDOC .NE. NULL) THEN C LOAD EXP. PROJ. ANGLES & ALIGNMENT PARAMETERS (ANGEXP) C FROM DOC. FILE (EXPANGDOC) CALL AP_GETANGA(ILIP,NIDI,0,EXPANGDOC,EXPPAT, & INPIC,LUNANG,8,ANGEXP,NGOTPAR,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ELSE ANGEXP = 0.0 ENDIF C INITIALIZE CCROT STATISTICS CCROTAVG = 0.0 IMPROVCCROT = 0 CCROTIMPROV = 0.0 IWORSECCROT = 0 CCROTWORSE = 0.0 ANGDIFAVG = 0.0 IBIGANGDIF = 0 C DISTRIBUTE EXPERIMENTAL IMAGES ALLOCATE(PSIZE(NPROCS), NBASE(NPROCS), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'DSGR_P: NBASE',2*NPROCS) GOTO 9999 ENDIF C FILL PSIZE WITH PARTITION LIMITS CALL SETPART(NIDI, PSIZE, NBASE) #ifdef MPI_DEBUG WRITE(6,111) NBASE(MYPID+1), MYPID 111 FORMAT(' DSGR_P: NBASE(MYPID+1): ', I5, ' MYPID: ', I5) CALL FLUSHFILE(6) CALL MPI_BARRIER(ICOMM,MPIERR) #endif NLOC = PSIZE(MYPID+1) ALLOCATE(ALOC(NSAM*NROW,NLOC), & ABUF(NSAM*NROW,PSIZE(1)), & PARTAB(15,NIDI), & PARTABLOC(15,NLOC), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = NSAM*NROW*NLOC + NSAM*NROW*PSIZE(1) + 15*(NIDI+NLOC) CALL ERRT(46,'DSGR_P: ALOC...',MWANT) RETURN ENDIF C ZERO THE ARRAYS ALOC = 0.0 ABUF = 0.0 PARTAB = 0.0 PARTABLOC = 0.0 #ifdef MPI_DEBUG TCOM0 = MPI_WTIME() #endif DO IPROC = 1, NPROCS NLOCP = PSIZE(IPROC) C === CALCULATE THE GLOBAL INDEX === IBEG = NBASE(IPROC) + 1 IEND = NBASE(IPROC) + NLOCP C ALTHOUGH THE FOLLOWING IS CALLED BY ALL PROCESSORS, C ONLY ONE PROCESSOR (MYPID=0) READS IMAGES INTO ABUF CALL AP_GETDATS(ILIP,NIDI,NSAM,NROW, & 1,EXPPAT,INPIC, IBEG,IEND, & ABUF, IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 CALL MPI_BARRIER(ICOMM, MPIERR) IF (IPROC .GT. 1) THEN IF (MYPID .EQ. 0) THEN #ifdef MPI_DEBUG WRITE(6,*)' DSGR_P: SENDING TO MYPID: ',IPROC-1 CALL FLUSHFILE(6) #endif CALL SEND_MPI('DSGR_P','ABUF', ABUF, NSAM*NROW*NLOCP, & 'R',IPROC-1,IPROC-1, ICOMM) ELSE IF (MYPID .EQ. IPROC-1) THEN C === SLAVES RECEIVE LOCAL PIECES === CALL RECV_MPI('DSGR_P','ALOC', ALOC, NSAM*NROW*NLOCP, & 'R', 0,0, ICOMM) #ifdef MPI_DEBUG WRITE(6,*) ' DSGR_P: RECEIVED BY MYPID;', MYPID CALL FLUSHFILE(6) #endif ENDIF ELSE IF (MYPID .EQ. 0) THEN C === SIMPLY COPY FROM ABUF TO ALOC === DO JLOC = 1, NLOCP DO IROW = 1, NSAM*NROW ALOC(IROW,JLOC) = ABUF(IROW,JLOC) ENDDO ENDDO ENDIF ENDDO ! END OF: DO IPROC = 1, NPROCS IF (ALLOCATED(ABUF)) DEALLOCATE(ABUF) #ifdef MPI_DEBUG_never TCOM1 = MPI_WTIME() IF (MYPID .EQ. 0) WRITE(6,440) TCOM1-TCOM0 440 FORMAT(' DSGR_P: DATA DIST TIME = ', 1PE11.3) #endif C CALCULATE DIMENSIONS FOR APRINGS CNS2 = NSAM / 2 + 1 CNR2 = NROW / 2 + 1 C THE FOLLOWING LOOP IS SIMULTANEOUSLY PERFORMED BY ALL PROCESSORS C ON DISTRIBUTED DATA c write(0,*)' dsgr_p; barrier 0,nloc,nprocs: ',nloc,nprocs,mypid CALL MPI_BARRIER(ICOMM,MPIERR) USEBCAST = .FALSE. DO JLOC = 1, NLOC ! --------------------------------------- IEXP = NBASE(MYPID+1) + JLOC C CONVERT EXP. ANGLE TO UNITARY DIRECTIONAL VECTORS (EXPDIR). CALL AP_GETSATA(ANGEXP(1,IEXP),EXPDIR,8,1,IRTFLG) C EXTRACT EXP. IMAGE POLAR COORD. RINGS, NORMALIZE & FFT THEM CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, ALOC(1,JLOC), & .FALSE.,MODE,NUMR,NRING,LCIRC, 0.0,FFTW_PLANS, & CIRCEXP,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C DETERMINE WHICH REF IMAGES ARE TO BE COMPARED NIMALCG = 0 IEND = NIMA CCROTD = -1.0D20 LIMITRANGE = RANGE .LT. 1.0 IF (LIMITRANGE) THEN C DETERMINE WHICH REF IMAGES ARE TO BE COMPARED DO IMI=1,NIMA C LOOP OVER ALL REF. IMAGES C ABS - DIRECTIONS AT 180 DEGREES ARE DIFFERENT C (- DO NOT CHECK MIRRORED) C DT NEAR 1.0 = NOT-MIRRORED, DT NEAR -1.0 = MIRRORED DT = (EXPDIR(1) * REFDIR(1,IMI) + & EXPDIR(2) * REFDIR(2,IMI) + & EXPDIR(3) * REFDIR(3,IMI)) DTABS = ABS(DT) IF (DTABS .GE. RANGE) THEN C NON-MIRRORED OR MIRRORED REF. IS WITHIN RANGE NIMALCG = NIMALCG + 1 LCG(NIMALCG) = IMI IF (DT .LT. 0) LCG(NIMALCG) = -IMI ENDIF ENDDO IF (NIMALCG .LE. 0) THEN C NO REF. IMAGE WITHIN COMPARISON ANGLE IREF = 0 GOTO 1000 ENDIF IEND = NIMALCG ENDIF C REF. IMAGES FOUND WITHIN COMPARISION RANGE IF (CIRCREF_IN_CORE) THEN C USE CIRCREF FOR REFERENCE RINGS DO IMIL=1,IEND IMI = IMIL IF (LIMITRANGE) IMI = ABS(LCG(IMIL)) IF ((CKMIRROR .AND. LIMITRANGE) .OR. & (.NOT. CKMIRROR)) THEN IF (.NOT. CKMIRROR) MIRRORED = .FALSE. IF (LIMITRANGE) MIRRORED = (LCG(IMIL) .LT. 0) C CHECK EITHER MIRRORED OR NON-MIRRORED POSITION CALL CROSRNG_EP_NEW(CIRCREF(1,IMI),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TT,MIRRORED,FFTW_PLANS(1)) ELSE C CHECK BOTH NON-MIRRORED & MIRRORED POSITIONS CALL CROSRNG_MSP_NEW(CIRCREF(1,IMI),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TOTMIR(IMI),ROTANGM(IMI), TT,FFTW_PLANS(1)) ENDIF ENDDO ! END OF: DO IMIL=1,IEND ELSE C USE REFERENCE RINGS FILE (MIGHT BE AN INCORE FILE) DO IMIL=1,IEND !LOOP OVER ALL REFERENCE IMAGES IMI = IMIL IF (LIMITRANGE) IMI = ABS(LCG(IMIL)) C FIND THREAD NUMBER ITHREAD = MOD((IMIL-1),NUMTH) + 1 C FILL CIRCREF FROM REFERENCE RINGS FILE CALL REDLIN(LUNRING,CIRCREF(1,ITHREAD),LCIRC,IMI) IF ((CKMIRROR .AND. LIMITRANGE) .OR. & (.NOT. CKMIRROR)) THEN IF (.NOT. CKMIRROR) MIRRORED = .FALSE. IF (LIMITRANGE) MIRRORED = (LCG(IMIL) .LT. 0) C CHECK EITHER MIRRORED OR NON-MIRRORED POSITION CALL CROSRNG_EP_NEW(CIRCREF(1,ITHREAD),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TT,MIRRORED,FFTW_PLANS(1)) ELSE C CHECK BOTH NON-MIRRORED & MIRRORED POSITIONS CALL CROSRNG_MSP_NEW(CIRCREF(1,ITHREAD),CIRCEXP, & LCIRC,NRING, MAXRIN,NUMR, & TOTMIN(IMI),ROTANGT(IMI), & TOTMIR(IMI),ROTANGM(IMI), & TT,FFTW_PLANS(1)) ENDIF ENDDO ! END OF: DO IMIL=1,IEND ENDIF C LOOP OVER ALL RELEVANT REF. IMAGES IREF = 0 CCROTD = -1.0D20 DO IMIL=1,IEND IMI = IMIL IF (NIMALCG .GT. 0) IMI = ABS(LCG(IMIL)) IF (TOTMIN(IMI) .GE. CCROTD) THEN C GOOD MATCH WITH TOTMIN (MIRRORED OR NOT) POSITION CCROTD = TOTMIN(IMI) RANGNEW = ROTANGT(IMI) MIRRORNEW = (LIMITRANGE .AND. (LCG(IMIL) .LT. 0)) IREF = IMI ENDIF IF (CKMIRROR .AND. .NOT. LIMITRANGE) THEN C HAVE TO COMPARE WITH MIRRORED POSITION IF (TOTMIR(IMI) .GE. CCROTD) THEN C GOOD MATCH, MIRRORED POSITION IS BETTER CCROTD = TOTMIR(IMI) RANGNEW = ROTANGM(IMI) MIRRORNEW = .TRUE. IREF = IMI ENDIF ENDIF ! END OF: IF (CKMIRROR) ENDDO ! END OF: DO IMIL=1,IEND 1000 RANGNEW = (RANGNEW-1) / MAXRIN * DIVAS CCROT = CCROTD IMGEXP = ILIP(IEXP) PEAKV = 0.0 XSHNEW = 0.0 YSHNEW = 0.0 IF (IREF .LE. 0) THEN C NO NEARBY REFERENCE IMAGE IMGREF = 0 C IREFT IS FOR REFDIR INDEX IREFT = 1 ELSE IMGREF = ILIST(IREF) C IREFT IS FOR REFDIR INDEX IREFT = IREF ENDIF NPROJ = NIMA IF (LIMITRANGE) NPROJ = NIMALCG C AP_END RETURNS PARTABLOC PARAMETERS FOR THIS IMAGE CALL AP_END(IEXP,IMGEXP,IMGREF, & ANGREF(1,IREFT),REFDIR(1,IREFT), & ANGEXP(1,IEXP),EXPDIR,ISHRANGE, & GOTREFANG, NGOTPAR,NSAM,NROW,CCROT,PEAKV, & RANGNEW,XSHNEW,YSHNEW,MIRRORNEW,EXPPAT,REFPAT, & NPROJ, CTYPE, ALOC(1,JLOC), LUNDOC, PARTABLOC(1,JLOC)) CCROTAVG = CCROTAVG + CCROT IF (NGOTPAR .GE. 8) THEN C COMPILE CCROT CHANGE STATISTICS ANGDIF = PARTABLOC(10,JLOC) IF (ANGDIF .GT. ANGDIFTHR)IBIGANGDIF = IBIGANGDIF + 1 CCROTLAS = ANGEXP(8,IEXP) ANGDIFAVG = ANGDIFAVG + PARTABLOC(10,JLOC) CCROTAVG = CCROTAVG + CCROT IF (CCROT .GE. CCROTLAS) THEN IMPROVCCROT = IMPROVCCROT + 1 CCROTIMPROV = CCROTIMPROV + CCROT ELSE IWORSECCROT = IWORSECCROT + 1 CCROTWORSE = CCROTWORSE + CCROT ENDIF ENDIF ENDDO ! END OF: DO JLOC = 1, NLOC ---------------------- c write(0,*) ' dsgr_p, barrier 2, mypid: ' ,mypid CALL MPI_BARRIER(ICOMM,MPIERR) USEBCAST = .TRUE. C COLLECT ALL ALIGNMENT PARAMETERS INTO PARTAB DO IPROC = 1, NPROCS PSIZE(IPROC) = 15 * PSIZE(IPROC) NBASE(IPROC) = 15 * NBASE(IPROC) ENDDO CALL MPI_ALLGATHERV(PARTABLOC, PSIZE(MYPID+1), & MPI_REAL, PARTAB, PSIZE, NBASE, & MPI_REAL, ICOMM, MPIERR) C WRITE IS SYNCHRONIZED WITHIN LUNDOCWRTDAT IF (LUNDOC .GT. 0) THEN C SAVE IN ALIGNMENT DOC FILE C <,<,<, MIR-REF#,IMG#,INPLANE<, SX,SY,NPROJ, C