C++********************************************************************* C C DSFS_P.F C CALL flush fixed Jul 2006 ArDean Leith C NUMREF BUG Mar 2008 ArDean Leith C FFTW_PLANS MAR 08 ARDEAN LEITH C APRINGS_ONE_NEW MAR 08 ARDEAN LEITH 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 DSFS_P(ILIST,NUMREF,IMGLST,NUMEXP, LSAM,LROW,NR,LENTT, C NRING,LCIRC,NUMR,CIRCREF,CIRCREF_IN_CORE, C MODE,SCRFILE,REFPAT,EXPPAT,FFTW3PLAN) C C PURPOSE: FIND ROTATIONAL 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_MS. CHECK BOTH STRAIGHT AND MIRRORED C POSITIONS BY CALCULATING X*CONJG(Y) AND CONJG(X)*CONJG(Y) C SINCE Y WERE PRE-WEIGHTED THE RESULTS IS ALREADY CORRECT. CC C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* SUBROUTINE DSFS_P(ILIST,NUMREF,IMGLST,NUMEXP, & LSAM,LROW,NR,LENTT, & NRING,LCIRC,NUMR,CIRCREF,CIRCREF_IN_CORE, & MODE,SCRFILE,REFPAT,EXPPAT,FFTW_PLANS) INCLUDE 'CMLIMIT.INC' INCLUDE 'CMBLOCK.INC' INTEGER,PARAMETER :: NLIST=7 INTEGER, DIMENSION(NUMREF) :: ILIST INTEGER, DIMENSION(NUMEXP) :: IMGLST INTEGER, DIMENSION(3,NRING) :: NUMR REAL, DIMENSION(LCIRC,NUMREF) :: CIRCREF CHARACTER (LEN=*) :: SCRFILE CHARACTER (LEN=*) :: REFPAT,EXPPAT CHARACTER(LEN=1) :: MODE REAL, ALLOCATABLE, DIMENSION(:,:) :: X CHARACTER(LEN=1) :: NULL LOGICAL :: CIRCREF_IN_CORE DOUBLE PRECISION :: EAV INTEGER *8 :: FFTW_PLANS(*) C AUTOMATIC ARRAYS REAL, DIMENSION(NLIST) :: DLIST DOUBLE PRECISION, DIMENSION(LENTT) :: TT REAL, DIMENSION(LCIRC) :: CIROLD, CIRC DOUBLE PRECISION, DIMENSION(NUMREF) :: TOTMIN,TOTMIR REAL, DIMENSION(NUMREF) :: TOT,TMT DATA INPIC/77/,NDOC/55/,NSCF/50/ #ifdef USE_MPI INCLUDE 'mpif.h' INTEGER IERR, COMM, MYPID, NPROCS, RC, TAG INTEGER MPISTAT(MPI_STATUS_SIZE) INTEGER NPRJLOC, NREM, JLOC, IPROC, NLOC, MASTER, & IROW, JCOL, IBEG INTEGER, ALLOCATABLE, DIMENSION(:) :: PSIZE INTEGER, ALLOCATABLE, DIMENSION(:) :: NBASE REAL, ALLOCATABLE, DIMENSION(:,:) :: DLISTGLB, DLISTLOC REAL, ALLOCATABLE, DIMENSION(:,:) :: CIRLOC REAL, ALLOCATABLE, DIMENSION(:,:) :: CIRBUF C COMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(COMM, MYPID, IERR) CALL MPI_COMM_SIZE(COMM, NPROCS, IERR) #else MYPID = -1 #endif NULL = CHAR(0) C ZERO DLIST ARRAY DLIST = 0.0 MAXRIN = NUMR(3,NRING) CALL APMASTER_1(MODE,DIVAS,NR,NUMTH,LSAM,LROW,NSAM,NROW, & TT,LENTT) ALLOCATE(X(NSAM,NROW),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = NSAM*NROW CALL ERRT(46,'X',MWANT) GOTO 9999 ENDIF #ifdef SP_MP IF (.NOT. CIRCREF_IN_CORE) THEN WRITE(NOUT,*) ' SETTING OMP THREADS: 2' CALL SETTHREADS(2) ENDIF #endif CALL REG_GET_USED(NSEL_USED) LQ = LROW/2+1 LR1 = (NROW-1)/2 LR2 = LQ+LR1 LR1 = LQ-LR1 LQ = LSAM/2+1 LS1 = (NSAM-1)/2 LS2 = LQ+LS1 LS1 = LQ-LS1 C READ REFERENCE IMAGES INTO REFERENCE RINGS (CIRCREF) ARRAY IF POSSIBLE CALL APRINGS_NEW(ILIST,NUMREF, LSAM,LROW, & NRING,LCIRC,NUMR,MODE,FFTW_PLANS, & REFPAT,INPIC,CIRCREF,CIRCREF_IN_CORE, & NSCF,SCRFILE,IRTFLG) C CALCULATE DIMENSIONS FOR APRINGS CNS2 = NSAM / 2 + 1 CNR2 = NROW / 2 + 1 #ifdef USE_MPI c --- START MPI VERSION ----------------------------------------- MASTER = 0 ALLOCATE(PSIZE(NPROCS), STAT=IRTFLG) IF ( IRTFLG. NE. 0) THEN WRITE(6,*) 'DSFS_P: FAILED TO ALLOCATE PSIZE' STOP ENDIF ALLOCATE(NBASE(NPROCS), STAT=IRTFLG) IF ( IRTFLG .NE. 0) THEN WRITE(6,*) 'DSFS_P: FAILED TO ALLOCATE NBASE' STOP ENDIF C C PARTITION THE PARTICLE IMAGES C THE FOURIER RINGS WILL BE STORED IN CIRLOC() C CALL SETPART(NUMEXP, PSIZE, NBASE) NPRJLOC = PSIZE(MYPID+1) ALLOCATE(CIRLOC(LCIRC,NPRJLOC), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN WRITE(6,*) 'DSFS_P: FAILED TO ALLOCATE CIRLOC' STOP ENDIF #ifdef MPI_DEBUG WRITE(6,111) PSIZE(MYPID+1), NBASE(MYPID+1), MYPID 111 FORMAT(1X, 'DSFS_P: PSIZE = ', I5, ' NBASE = ', I5, & ' MYPID = ', I3) #endif C C === ALLOCATE A BUFFER TO READ IMAGE === C ALLOCATE(CIRBUF(LCIRC, PSIZE(1)), STAT=IRTFLG) IF ( IRTFLG .NE. 0 ) THEN WRITE(6,*) 'DSFS_P: FAILED TO ALLOCATE CIRBUF' STOP ENDIF DO IPROC = 1, NPROCS NLOC = PSIZE(IPROC) C ALTHOUGH AP_GETDAT1P IS CALLED BY ALL PROCESSORS, C ONLY ONE PROCESSOR (MYPID=0) READS IMAGES INTO CIRBUF DO JLOC = 1, NLOC ITI = NBASE(IPROC) + JLOC c LOAD WINDOW FROM SAMPLE IMAGES INTO ARRAY X CALL AP_GETDAT1P(IMGLST,NUMEXP,LSAM,LROW,NSAM,NROW, & NUMTH,EXPPAT,INPIC, ITI,ITI, & LR1,LR2,LS1,LS2, X, & IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C ONLY THE MASTER NEEDS TO DO THE FOLLOWING IF (MYPID .LE. 0) THEN C EXTRACT EXP. IMAGE POLAR COORD. RINGS, NORMALIZE & C FFT THEM, DANGER THIS CALLS FRNGS (NOT OMP ||) CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, X, .FALSE., & MODE,NUMR,NRING,LCIRC, 0.0,FFTW_PLANS, & CIRBUF(1,JLOC),IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ENDIF ENDDO C SEND THE FOURIER RINGS OFF TO OTHER PROCESSORS IF (IPROC .GT. 1) THEN IF (MYPID .LE. 0) THEN #ifdef MPI_DEBUG WRITE(6,112) IPROC-1 CALL FLUSHFILE(6) 112 FORMAT(' DSFS_P: SENDING TO PID = ', I3) #endif CALL MPI_SEND(CIRBUF , NLOC*LCIRC, MPI_REAL, & IPROC-1, IPROC-1 , COMM , & IERR) ELSE IF (MYPID .EQ. IPROC-1) THEN C === SLAVES RECEIVE LOCAL PIECES === CALL MPI_RECV(CIRLOC, NPRJLOC*LCIRC, MPI_REAL, & 0 , MPI_ANY_TAG , COMM , & MPISTAT , IERR) IF (IERR .NE. 0) THEN WRITE(6,*) ' RECV FAILED' STOP ENDIF #ifdef MPI_DEBUG WRITE(6,113) MYPID CALL FLUSHFILE(6) 113 FORMAT(' DSFS_P: MYPID = ', I3, ' RECEIVED') #endif ENDIF ELSE IF (MYPID .EQ. 0) THEN C === SIMPLY COPY FROM CIRBUF TO CIRLOC === DO JCOL = 1, NLOC DO IROW = 1, LCIRC CIRLOC(IROW,JCOL) = CIRBUF(IROW,JCOL) END DO END DO END IF END DO CALL MPI_BARRIER(COMM, IERR) DEALLOCATE(CIRBUF) C ALLOCATE(DLISTGLB(7,NUMEXP), DLISTLOC(7,NPRJLOC), & STAT = IRTFLG) IF (IRTFLG .NE. 0) THEN WRITE(6,*) 'DSFS_P: FAILED TO ALLOCATE DLISTGLB, DLISTLOC' STOP endif DLISTGLB = 0.0 DLISTLOC = 0.0 C === PERFORMED SIMULTANEOUSLY BY ALL PROCESSORS === DO ITI = 1, NPRJLOC EAV = -1.0D20 IF (CIRCREF_IN_CORE) THEN DO IMI=1,NUMREF CALL CROSRNG_MSP_NEW(CIRCREF(1,IMI) , CIRLOC(1,ITI), & LCIRC , NRING , & MAXRIN , NUMR , & TOTMIN(IMI), TOT(IMI) , & TOTMIR(IMI), TMT(IMI) , & TT,FFTW_PLANS(1)) ENDDO ELSE C === REFERENCE RINGS OUT OF CORE === IF (MYPID .LE. 0) REWIND NSCF DO IMI=1,NUMREF IF (MYPID .LE. 0) THEN READ(NSCF) CIRC ENDIF CALL MPI_BCAST(CIRC,LCIRC,MPI_REAL,0,COMM,IERR) CALL CROSRNG_MSP_NEW(CIRC,CIRLOC(1,ITI),LCIRC,NRING, & MAXRIN,NUMR,TOTMIN(IMI),TOT(IMI),TOTMIR(IMI), & TMT(IMI),TT,FFTW_PLANS(1)) ENDDO ENDIF C === PULL OUT THE RESULTS === DO IMI=1,NUMREF IF ( TOTMIN(IMI) .GE. EAV .AND. & TOTMIN(IMI) .GT. TOTMIR(IMI) ) THEN EAV = TOTMIN(IMI) IDI = ILIST(IMI) RANG = TOT(IMI) ELSEIF ( TOTMIR(IMI) .GE. EAV ) THEN EAV = TOTMIR(IMI) IDI = -ILIST(IMI) RANG = TMT(IMI) ENDIF ENDDO IGLB = ITI + NBASE(MYPID + 1) DLISTLOC(1,ITI) = IGLB DLISTLOC(2,ITI) = IDI DLISTLOC(3,ITI) = EAV RANG = (RANG-1)/MAXRIN*DIVAS DLISTLOC(4,ITI) = RANG C DLIST 5&6 ARE PERMANENTLY SET TO ZERO (THIS IS TO KEEP C SAME FORMAT AS 'AP MQ' COMMAND) DLISTLOC(7,ITI) = IMGLST(IGLB) ENDDO CALL MPI_BARRIER(COMM,IERR) C === MERGE ALL DLISTLOC BY A GLOBAL SUM === DO IPROC = 1, NPROCS PSIZE(IPROC) = 7*PSIZE(IPROC) NBASE(IPROC) = 7*NBASE(IPROC) ENDDO CALL MPI_ALLGATHERV(DLISTLOC, PSIZE(MYPID+1), & MPI_REAL, DLISTGLB, PSIZE, NBASE, & MPI_REAL, COMM , IERR) DO ITI = 1, NUMEXP CALL LUNDOCWRTDAT(NDOC,ITI,DLISTGLB(2,ITI),NLIST-1, & IRTFLG) ENDDO IF (ALLOCATED(DLISTGLB)) DEALLOCATE(DLISTGLB) IF (ALLOCATED(DLISTLOC)) DEALLOCATE(DLISTLOC) IF (ALLOCATED(PSIZE)) DEALLOCATE(PSIZE) IF (ALLOCATED(NBASE)) DEALLOCATE(NBASE) IF (ALLOCATED(NBASE)) DEALLOCATE(NBASE) IF (ALLOCATED(CIRLOC)) DEALLOCATE(CIRLOC) CALL MPI_BARRIER(COMM,IERR) CALL FLUSHRESULTS c --- END OF MPI VERSION ---------------------------------- #else DO ITI=1,NUMEXP c write(6,*) ' --scrfile:',scrfile c write(6,*) ' --refpat:',refpat c write(6,*) ' --exppat:',exppat c write(6,*) ' --NUMEXP:',NUMEXP c write(6,*) ' --iti:',iti C LOAD WINDOW FROM SAMPLE IMAGES INTO ARRAY X CALL AP_GETDAT(IMGLST,NUMEXP,LSAM,LROW,NSAM,NROW, & NUMTH,EXPPAT,INPIC, ITI,ITI, & LR1,LR2,LS1,LS2, X, IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C EXTRACT EXP. IMAGE POLAR COORD. RINGS, NORMALIZE & C FFT THEM, DANGER THIS CALLS FRNGS , OMP || FRNG CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, X, .TRUE., & MODE,NUMR,NRING,LCIRC, 0.0,FFTW_PLANS, & CIROLD,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 EAV = -1.0D20 IF (CIRCREF_IN_CORE) THEN c$omp parallel do private(imi) DO IMI=1,NUMREF CALL CROSRNG_MSP_NEW(CIRCREF(1,IMI),CIROLD, & LCIRC,NRING, & MAXRIN,NUMR,TOTMIN(IMI),TOT(IMI),TOTMIR(IMI), & TMT(IMI),TT,FFTW_PLANS(1)) ENDDO ELSE REWIND NSCF DO IMI=1,NUMREF READ(NSCF) CIRC CALL CROSRNG_MSP_NEW(CIRC,CIROLD, & LCIRC,NRING, & MAXRIN,NUMR,TOTMIN(IMI),TOT(IMI),TOTMIR(IMI), & TMT(IMI),TT,FFTW_PLANS(1)) ENDDO ENDIF DO IMI=1,NUMREF IF (TOTMIN(IMI) .GE. EAV .AND. & TOTMIN(IMI) .GT. TOTMIR(IMI)) THEN EAV = TOTMIN(IMI) IDI = ILIST(IMI) RANG = TOT(IMI) ELSEIF (TOTMIR(IMI) .GE. EAV) THEN EAV = TOTMIR(IMI) IDI = -ILIST(IMI) RANG = TMT(IMI) ENDIF ENDDO DLIST(1) = ITI DLIST(2) = IDI DLIST(3) = EAV RANG = (RANG-1 )/ MAXRIN*DIVAS DLIST(4) = RANG C DLIST 5&6 ARE SET TO ZERO (KEEP SAME FORMAT AS 'AP MQ') DLIST(7) = IMGLST(ITI) C CALL SAVDN1(NDOC,CDUM,DLIST,NLIST, 1,0) CALL LUNDOCWRTDAT(NDOC,ITI,DLIST(2),NLIST-1,IRTFLG) CALL FLUSHFILE(NDOC) ENDDO #endif 9999 IF (.NOT. CIRCREF_IN_CORE) CLOSE(NSCF) IF (MYPID .LE. 0) CLOSE(NDOC) #ifdef USE_MPI call MPI_BARRIER(comm,ierr) #endif IF (ALLOCATED(X)) DEALLOCATE(X) END