C++********************************************************************* C C WIW32D.F 01/10/00 C ANGLES IN HEADER JULY 01 ArDean Leith C OPFILE 'U' FEB 02 ArDean Leith C BESSELS 05/03/02 P.A.Penczek C OPFILEC FEB 03 ARDEAN LEITH C BUILDM PARAMETERS JUL 03 ARDEAN LEITH C BUILDM PARAMETERS SEP 03 ARDEAN LEITH C MPI OCT 03 CHAO YANG C WIW32D_DL.F FEB 07 ARDEAN LEITH C FMRS_PLAN MAY 08 ARDEAN LEITH C C=********************************************************************** C=* From: SPIDER - MODULAR IMAGE PROCESSING SYSTEM * C=* Copyright (C)2002 & 2008, P. A. Penczek * C=* * C=* University of Texas - Houston Medical School * C=* * C=* Email: pawel.a.penczek@uth.tmc.edu * 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 WIW32D C C PURPOSE: FOR 'BP 32F' BACK PROJECTION - 3D, ODD-EVEN, INTERPOLATED C IN FOURIER SPACE ||* C AND FOR: 'BPD 32F' IMPROVED WITH LESS MEMORY REQUIRED AND C MORE STACK USAGE. C23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 C--********************************************************************* SUBROUTINE WIW32D INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' INCLUDE 'F90ALLOC.INC' C DOC FILE POINTERS REAL, DIMENSION(:,:), POINTER :: ANGBUF, ANGSYM REAL, DIMENSION(:,:), ALLOCATABLE :: DM,SM COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: XE,XO REAL, ALLOCATABLE, DIMENSION(:,:,:) :: WE,WO REAL, DIMENSION(:), ALLOCATABLE :: TEMP LOGICAL :: ANGINDOC CHARACTER(LEN=1) :: NULL CHARACTER*80 FINPIC,FINFO,ANGDOC,FINPAT COMMON /F_SPEC/ FINPAT,NLET,FINPIC DATA IOPIC/98/,INPIC/99/ #ifdef USE_MPI INCLUDE 'mpif.h' INTEGER MYPID, COMM, MPIERR COMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(COMM, MYPID, MPIERR) #else MYPID = -1 #endif NULL = CHAR(0) NILMAX = NIMAX CALL FILELIST(.TRUE.,INPIC,FINPAT,NLET,INUMBR,NILMAX,NANG, & 'ENTER TEMPLATE FOR 2-D IMAGES',IRTFLG) IF (IRTFLG .NE. 0) RETURN CLOSE(INPIC) MAXNUM = MAXVAL(INUMBR(1:NANG)) C NANG - TOTAL NUMBER OF IMAGES IF (MYPID .LE. 0) WRITE(NOUT,2001) NANG 2001 FORMAT(' NUMBER OF IMAGES: ',I7) C RETRIEVE ARRAY WITH ANGLES DATA IN IT ANGINDOC = .TRUE. MAXXT = 4 MAXYT = MAXNUM CALL GETDOCDAT('ANGLES DOC',.TRUE.,ANGDOC,77,.FALSE.,MAXXT, & MAXYT,ANGBUF,IRTFLG) IF (IRTFLG .NE. 0) ANGINDOC = .FALSE. C RETRIEVE ARRAY WITH SYMMETRIES DATA IN IT MAXXS = 0 MAXSYM = 0 CALL GETDOCDAT('SYMMETRIES DOC',.TRUE.,ANGDOC,77,.TRUE.,MAXXS, & MAXSYM,ANGSYM,IRTFLG) IF (IRTFLG .NE. 0) MAXSYM=1 C OPEN FIRST IMAGE FILE TO DETERMINE NSAM, NROW, NSL CALL FILGET(FINPAT,FINPIC,NLET,INUMBR(1),IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINPIC,INPIC,'O',IFORM,NSAM,NROW,NSL, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CLOSE(INPIC) N2 = 2*NSAM LSD = N2+2-MOD(N2,2) NMAT = LSD*N2*N2 IF (ANGINDOC) THEN ALLOCATE(DM(9,NANG), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN MEMWANT = 9 * NANG CALL ERRT(46,'BP 32F, DM',MEMWANT) GOTO 999 ENDIF CALL BUILDM(INUMBR,DM,NANG,ANGBUF(1,1),.FALSE.,SSDUM, & .FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 IF (ASSOCIATED(ANGBUF)) DEALLOCATE(ANGBUF) ELSE ALLOCATE(DM(9,1), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'BP 32F, DM',IER) GOTO 999 ENDIF ENDIF IF (MAXSYM .GT. 1) THEN ALLOCATE(SM(9,MAXSYM), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN MEMWANT = 9 * MAXSYM CALL ERRT(46,'BP 3F, SM',MEMWANT) GOTO 999 ENDIF CALL BUILDS(SM,MAXSYM,ANGSYM(1,1),IRTFLG) IF (ASSOCIATED(ANGSYM)) DEALLOCATE(ANGSYM) ELSE ALLOCATE(SM(1,1), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'BP 3F, SM(1,1)',IER) GOTO 999 ENDIF ENDIF ALLOCATE(XE(0:N2/2,N2,N2), & WE(0:N2/2,N2,N2), & WO(0:N2/2,N2,N2), & XO(0:N2/2,N2,N2), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN C SOME ARRAYS ARE COMPLEX SO 6 NOT 4 MEMWANT = 6 * ((N2/2)+1) * N2 * N2 CALL ERRT(46,'BP 32F; XE, WE, WO, & XO',MEMWANT) GOTO 999 ENDIF C CREATE FFTW3 PLAN FOR 3D FFT ON XE USING ALL THREADS CALL FMRS_PLAN(.TRUE.,XE,N2,N2,N2, 0,-1,IRTFLG) IF (IRTFLG .NE. 0) RETURN CALL WIW32DQ(NSAM,XE,WE,XO,WO, & LSD,N2,N2/2,INUMBR,DM,NANG,SM,MAXSYM,ANGINDOC,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL FILERD(FINPIC,NLETI,NULL,'RECONSTRUCTED 3D',IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL FILERD(FINPAT,NLETI,NULL,'RECONSTRUCTED EVEN 3D',IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL FILERD(FINFO,NLETI,NULL,'RECONSTRUCTED ODD 3D',IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 C STORE FOURIER XE TEMPORARILY IFORM = 3 MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINPIC,IOPIC,'N',IFORM,LSD,N2,N2, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL WRITEV(IOPIC,XE,LSD,N2,LSD,N2,N2) IF (MYPID .LE. 0) CLOSE(IOPIC) C CBE AND CBO (XE & XO) ARE LSD x N2 x N2 C CWE AND CWO (WE & WO) ARE (LSD/2) x N2 x N2 LSDD2 = LSD / 2 IFORM = 3 MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINFO,IOPIC,'N',IFORM,LSDD2,N2,N2, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 C STORE WEIGHT WE TEMPORARILY CALL WRITEV(IOPIC,WE,LSD/2,N2,LSD/2,N2,N2) IF (MYPID .LE. 0) CLOSE(IOPIC) CALL NRMW2(XE,WE,N2/2,N2) CALL WINDKB2(XE,XE,NSAM,LSD,N2) C NOW XE IS READY, SYMMETRIZE IF NECESSARY C ADDITIONAL SYMMETRIZATION OF THE VOLUME XE IN REAL SPACE 05/03/02 IF (MAXSYM .GT. 1) THEN ALLOCATE (TEMP(NSAM*NSAM*NSAM), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN MEMWANT = NSAM * NSAM * NSAM CALL ERRT(46,'BP 32F, TEMP',MEMWANT) GOTO 999 ENDIF CALL COP(XE,TEMP,NSAM*NSAM*NSAM) c$omp parallel do private(i,j,k) DO K=1,N2 DO J=1,N2 DO I=0,N2/2 XE(I,J,K) = CMPLX(0.0,0.0) ENDDO ENDDO ENDDO IF (MOD(NSAM,2) .EQ. 0) THEN KNX = NSAM/2-1 ELSE KNX = NSAM/2 ENDIF KLX = -NSAM/2 CALL SYMVOL(TEMP,XE,KLX,KNX,KLX,KNX,KLX,KNX,SM,MAXSYM) ENDIF IFORM = 3 MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINPAT,IOPIC,'U',IFORM,NSAM,NSAM,NSAM, & MAXIM,' ',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL WRITEV(IOPIC,XE,NSAM,NSAM,NSAM,NSAM,NSAM) IF (MYPID .LE. 0) CLOSE(IOPIC) IFORM = 3 MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINPIC,IOPIC,'O',IFORM,LSD,N2,N2, & MAXIM,' ',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL READV(IOPIC,XE,LSD,N2,LSD,N2,N2) CLOSE(IOPIC) IFORM = 3 MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINFO,IOPIC,'O',IFORM,LSDD2,N2,N2, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL READV(IOPIC,WE,LSD/2,N2,LSD/2,N2,N2) CLOSE(IOPIC) C Add E+O CALL ADDADA(XE,XO,NMAT) CALL ADDADA(WE,WO,NMAT/2) CALL NRMW2(XE,WE,N2/2,N2) CALL WINDKB2(XE,XE,NSAM,LSD,N2) C TOTAL VOLUME IS READY C ADDITIONAL SYMMETRIZATION OF THE VOLUME TOTAL IN REAL SPACE 05/03/02 IF (MAXSYM .GT. 1) THEN CALL COP(XE,TEMP,NSAM*NSAM*NSAM) c$omp parallel do private(i,j,k) DO K=1,N2 DO J=1,N2 DO I=0,N2/2 XE(I,J,K)=CMPLX(0.0,0.0) ENDDO ENDDO ENDDO CALL SYMVOL(TEMP,XE,KLX,KNX,KLX,KNX,KLX,KNX,SM,MAXSYM) ENDIF MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINPIC,IOPIC,'U',IFORM,NSAM,NSAM,NSAM, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL WRITEV(IOPIC,XE,NSAM,NSAM,NSAM,NSAM,NSAM) IF (MYPID .LE. 0) CLOSE(IOPIC) CALL NRMW2(XO,WO,N2/2,N2) CALL WINDKB2(XO,XO,NSAM,LSD,N2) C XO VOLUME IS READY C ADDITIONAL SYMMETRIZATION OF THE VOLUME XO IN REAL SPACE 05/03/02 IF (MAXSYM .GT. 1) THEN CALL COP(XO,TEMP,NSAM*NSAM*NSAM) c$omp parallel do private(i,j,k) DO K=1,N2 DO J=1,N2 DO I=0,N2/2 XO(I,J,K) = CMPLX(0.0,0.0) ENDDO ENDDO ENDDO CALL SYMVOL(TEMP,XO,KLX,KNX,KLX,KNX,KLX,KNX,SM,MAXSYM) ENDIF MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINFO,IOPIC,'U',IFORM,NSAM,NSAM,NSAM, & MAXIM,' ',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL WRITEV(IOPIC,XO,NSAM,NSAM,NSAM,NSAM,NSAM) IF (MYPID .LE. 0) CLOSE(IOPIC) 999 CONTINUE IF (ALLOCATED(DM)) DEALLOCATE(DM) IF (ALLOCATED(SM)) DEALLOCATE(SM) IF (ALLOCATED(XE)) DEALLOCATE(XE) IF (ALLOCATED(WE)) DEALLOCATE(WE) IF (ALLOCATED(XO)) DEALLOCATE(XO) IF (ALLOCATED(WO)) DEALLOCATE(WO) IF (ALLOCATED(TEMP)) DEALLOCATE(TEMP) IF (ASSOCIATED(ANGBUF)) DEALLOCATE(ANGBUF) IF (ASSOCIATED(ANGSYM)) DEALLOCATE(ANGSYM) END C ---------------- WIW32DQ ------------------------------------- SUBROUTINE WIW32DQ(NS,XE,WE,XO,WO,LSD,N,N2, & INUMBRT,DM,NANG,SM,MAXSYM,ANGINDOC,IRTFLG) C NOTE: STUPID TRANSFORM OF N2-->N AND N2/2-->N2 !!!!al INCLUDE 'CMLIMIT.INC' LOGICAL :: ANGINDOC LOGICAL, DIMENSION(:), ALLOCATABLE :: RANDLIST REAL, DIMENSION(:,:), ALLOCATABLE :: PROJ COMPLEX, DIMENSION(:,:), ALLOCATABLE :: BI REAL, DIMENSION(4) :: ANGBUF DIMENSION WE(0:N2,N,N),WO(0:N2,N,N) COMPLEX XE(0:N2,N,N),XO(0:N2,N,N) DIMENSION INUMBRT(NANG) DIMENSION DM(3,3,NANG),SM(3,3,MAXSYM),DMS(3,3) CHARACTER*80 FINPIC,FINPAT COMMON /F_SPEC/ FINPAT,NLET,FINPIC DOUBLE PRECISION PI LOGICAL ITMP PARAMETER (LTAB=4999) COMMON /TABS/ LN2,FLTB,TABI(0:LTAB) C IN THIS VERSION THE ORDER OF THE BESSEL FUNCTION IS MMM=1 COMMON /BESSEL_PARAM/ ALPHA,AAAA,NNN DATA INPROJ/99/ PARAMETER (QUADPI = 3.141592653589793238462643383279502884197) PARAMETER (TWOPI = 2*QUADPI) #ifdef USE_MPI c This MPI version is memory intensive. c It requires 6 copies of the 3-D volume, c 2-D images are read into memory and distributed. c Each processor will hold roughly nang/nproc 2-D images. c Memory requirement may be reduced in a future version. INCLUDE 'mpif.h' INTEGER MYPID, COMM, MPIERR, NPROCS, K3 INTEGER ISTAT(MPI_STATUS_SIZE) INTEGER, ALLOCATABLE, DIMENSION(:) :: PSIZE, NBASE INTEGER NANGLOC, NREM, IPROC, ISAM, JROW, JGLB, JLOC REAL , ALLOCATABLE, DIMENSION(:,:,:) :: PRJLOC REAL , ALLOCATABLE, DIMENSION(:,:,:) :: PRJBUF REAL , ALLOCATABLE, DIMENSION(:,:,:) :: WELOC, WOLOC COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: XELOC, XOLOC COMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(COMM, MYPID, MPIERR) CALL MPI_COMM_SIZE(COMM, NPROCS, MPIERR) #else MYPID = -1 #endif C K=6 LN=5 LN2=LN/2 C GENERALIZED KAISER-BESSEL WINDOW ACCORDING TO LEWITT C M=NS, N=N R=NS/2 V=REAL(LN-1)/2.0/REAL(N) ALPHA=6.5 C AAAA=0.0079 C AAAA=1.1*V AAAA=0.9*V NNN=3 C mmm=1 C GENERATE TABLE WITH INTERPOLANTS C B0=(SQRT(ALPHA)**mmm)*BESSI(mmm,ALPHA) B0=SQRT(ALPHA)*BESI1(ALPHA) C if(abs(b0-tempo).gt.1.0e-5) print *,'B0' FLTB=REAL(LTAB)/REAL(LN2+1) C CANNOT BE PARALLEL AS THERE ARE DATA STATEMENTS IN BESSI cc$omp parallel do private(i,s,x),shared(mmm) DO I=0,LTAB S=REAL(I)/FLTB/N IF (S.LE.AAAA) THEN X=SQRT(1.0-(S/AAAA)**2) TABI(I)= SQRT(ALPHA*X)*BESI1(ALPHA*X)/B0 ELSE TABI(I)=0.0 ENDIF ENDDO c$omp parallel do private(i,j,k) DO K=1,N DO J=1,N DO I=0,N2 XE(I,J,K) = CMPLX(0.0,0.0) WE(I,J,K) = 0.0 XO(I,J,K) = CMPLX(0.0,0.0) WO(I,J,K) = 0.0 ENDDO ENDDO ENDDO ALLOCATE (PROJ(NS,NS), & BI(0:N2,N), & RANDLIST(NANG), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN MEMWANT = NS*NS + (N2 + 1) * N + NANG CALL ERRT(46,'BP 32F, PROJ, ...',MEMWANT) RETURN ENDIF C CREATE FFTW3 PLAN FOR 2D FFT ON BI USING ALL THREADS CALL FMRS_PLAN(.TRUE.,BI,N,N,1, 0,+1,IRTFLG) IF (IRTFLG .NE. 0) RETURN RANDLIST(1:NANG/2) = .TRUE. RANDLIST(NANG/2+1:NANG) = .FALSE. DO K=1,NANG CALL RANDOM_NUMBER(HARVEST=X) IORD = MIN(NANG,MAX(1,INT(X*NANG+0.5))) ITMP = RANDLIST(IORD) RANDLIST(IORD) = RANDLIST(K) RANDLIST(K) = ITMP ENDDO #ifdef SP_MP LN1 = LN+1 #endif #ifdef USE_MPI C --- BEGIN MPI VERSION --- C DISTRIBUTE PARTICLES TO PROCESSORS. C NANGLOC IS THE NUMBER OF PARTICLES ASSIGNED TO EACH PROCESSOR. ALLOCATE(PSIZE(NPROCS) ,NBASE(NPROCS), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MEMWANT = 2 * NPROCS CALL ERRT(46,'BP 3F, PSIZE & NBASE',MEMWANT) RETURN ENDIF C SETPART RETURNS THE SIZE OF THE LOCAL PIECE AND GLOBAL OFFSET. CALL SETPART(NANG, PSIZE, NBASE) NANGLOC = PSIZE(MYPID+1) C 2-D IMAGES ARE DISTRIBUTED AND HELD IN PRJLOC ON EACH PROCESSOR ALLOCATE(PRJBUF(NS,NS,PSIZE(1)), & PRJLOC(NS,NS,NANGLOC), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MEMWANT = NS*NS*PSIZE(1) + NS*NS*NANGLOC CALL ERRT(46,'BP 3F, PRJBUF, PRJLOC',MEMWANT) RETURN ENDIF C PROCESSOR 0 READS IMAGE FILES AND DISTRIBUTE THEM. C (THIS VERSION ASSUMES THAT THERE IS SUFFICIENT C MEMORY TO HOLD NANG/NPROCS IMAGES) DO IPROC = 1, NPROCS NANGLOC = PSIZE(IPROC) C READ IMAGES INTO THE BUFFER FIRST, THEN DISTRIBUTED C READ A BUNCH OF 2-D IMAGES DO JLOC = 1, NANGLOC JGLB = NBASE(IPROC) + JLOC CALL FILGET(FINPAT,FINPIC,NLET,INUMBR(JGLB),IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 MAXIM = 0 CALL OPFILEC(0 ,.FALSE., FINPIC, INPROJ, 'O' , & IFORM , NSAM , NSAM , NSL , MAXIM, & 'DUMMY',.FALSE., IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 CALL READV1P(INPROJ,PRJBUF(1,1,JLOC),NS,NS,NS,NS,1) CLOSE(INPROJ) ENDDO C IF (IPROC .GT. 1) THEN IF (MYPID .EQ. 0) THEN C SEND TO ANOTHER PROCESSOR #ifdef MPI_DEBUG WRITE(6,*) 'WIW32DQ: SENDING TO PID = ', IPROC-1 #endif CALL MPI_SEND(PRJBUF , NS*NS*NANGLOC, MPI_REAL, & IPROC-1, IPROC-1 , COMM , & MPIERR) IF (MPIERR .NE. 0) THEN WRITE(6,*) 'WIW32DQ: SEND ERROR!' STOP ENDIF ELSE IF (MYPID .EQ. IPROC-1) THEN C RECEIVE PROJECTION IMAGES FROM PROCESSOR 0 CALL MPI_RECV(PRJLOC, NS*NS*NANGLOC, MPI_REAL, & 0 , MPI_ANY_TAG , COMM , & ISTAT , MPIERR) IF (MPIERR .NE. 0) THEN WRITE(6,*) 'WIW32DQ: RECV FAILED' STOP ENDIF #ifdef MPI_DEBUG WRITE(6,111) MYPID 111 FORMAT(1X,'WIW32DQ: MYPID = ', I3, ' RECEIVED!') #endif ENDIF ELSE IF (MYPID .EQ. 0) THEN C KEEP FOR MYSELF DO JLOC = 1, NANGLOC DO ISAM = 1, NS DO JROW = 1, NS PRJLOC(ISAM,JROW,JLOC) & = PRJBUF(ISAM,JROW,JLOC) ENDDO ENDDO ENDDO ENDIF ENDDO IF (ALLOCATED(PRJBUF)) DEALLOCATE(PRJBUF) IF (.NOT. ANGINDOC) THEN C GET ANGLES FROM HEADER ANGBUF(1) = INUMBRT(K) CALL LUNGETVALS(INPROJ,IAPLOC+1,3,ANGBUF(2),IRTFLG) CALL BUILDM(INUMBRT,DM,1,ANGBUF,.FALSE.,SSDUM, & .FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ENDIF CLOSE(INPROJ) C PERFORM CALCULATIONS IN PARALLEL NOW ALLOCATE (WELOC(0:N2,N,N), XELOC(0:N2,N,N), & WOLOC(0:N2,N,N), XOLOC(0:N2,N,N), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MEMWANT = 4*(N2+1)*N*N CALL ERRT(46,'BP 3F, WELOC, XELOC, ...',MEMWANT) RETURN ENDIF DO K=1,N DO J=1,N DO I=0,N2 XELOC(I,J,K) = CMPLX(0.0,0.0) XOLOC(I,J,K) = CMPLX(0.0,0.0) WELOC(I,J,K) = 0.0 WOLOC(I,J,K) = 0.0 ENDDO ENDDO ENDDO C NANGLOC = PSIZE(MYPID+1) DO JLOC = 1, NANGLOC JGLB = NBASE(MYPID+1) + JLOC C PAD: PRJLOC TO SIZE: N CALL PADD2(PRJLOC(1,1,JLOC),NS,BI,LSD,N) C FOURIER TRANSFORM OF: BI INV = +1 CALL FMRS_2(BI,N,N,INV) DO J=1,N DO I=0,N2 BI(I,J)=BI(I,J)*(-1)**(I+J+1) ENDDO ENDDO DO ISYM=1,MAXSYM IF (MAXSYM .GT. 1) THEN C SYMMETRIES, MULTIPLY MATRICES DMS = MATMUL(SM(:,:,ISYM),DM(:,:,JGLB)) ELSE DMS = DM(:,:,JGLB) ENDIF IF (RANDLIST(K)) THEN DO J=-N2+1,N2 CALL ONELINE(J,N,N2,XEloc,WEloc,BI,DMS) ENDDO ELSE DO J=-N2+1,N2 CALL ONELINE(J,N,N2,XOloc,WOloc,BI,DMS) ENDDO ENDIF ENDDO ENDDO IF (ALLOCATED(PRJLOC)) DEALLOCATE(PRJLOC) C SUM UP X AND W FROM THE LOCAL PIECES (XLOC, WLOC) C RESIDING ON EACH PROCESSOR DO K3 = 1, N CALL MPI_ALLREDUCE(XELOC(1,1,K3), XE(1,1,K3), (N2+1)*N, & MPI_COMPLEX , MPI_SUM , COMM , & MPIERR) IF ( MPIERR. NE. 0 ) THEN WRITE(6,*) 'WIW32DQ: FAILED TO ALLREDUCE' STOP ENDIF CALL MPI_ALLREDUCE(XOLOC(1,1,K3), XO(1,1,K3), (N2+1)*N, & MPI_COMPLEX , MPI_SUM , COMM , & MPIERR) IF ( MPIERR. NE. 0 ) THEN WRITE(6,*) 'WIW32DQ: FAILED TO ALLREDUCE' STOP ENDIF CALL MPI_ALLREDUCE(WELOC(1,1,K3), WE(1,1,K3), (N2+1)*N, & MPI_REAL , MPI_SUM , COMM , & MPIERR) IF ( MPIERR. NE. 0 ) THEN WRITE(6,*) 'WIW32DQ: FAILED TO ALLREDUCE' STOP ENDIF CALL MPI_ALLREDUCE(WOLOC(1,1,K3), WO(1,1,K3), (N2+1)*N, & MPI_REAL , MPI_SUM , COMM , & MPIERR) IF ( MPIERR. NE. 0 ) THEN WRITE(6,*) 'WIW32DQ: FAILED TO ALLREDUCE' STOP ENDIF ENDDO IF (ALLOCATED(XOLOC)) DEALLOCATE (XOLOC) IF (ALLOCATED(XELOC)) DEALLOCATE (XELOC) IF (ALLOCATED(WOLOC)) DEALLOCATE (WOLOC) IF (ALLOCATED(WEloc)) DEALLOCATE (WEloc) #ifdef MPI_DEBUG WRITE(6,*) 'WIW32DQ: COMPLETED GLOBAL SUM, MYPID = ', MYPID #endif C --- END OF MPI VERSION --- #else DO K=1,NANG C print *,' Projection #',K C OPEN DESIRED FILE CALL FILGET(FINPAT,FINPIC,NLET,INUMBRT(K),IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINPIC,INPROJ,'O',IFORM,NSAM,NROW,NSL, & MAXIM,' ',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 CALL READV(INPROJ,PROJ,NS,NS,NS,NS,1) IF (.NOT. ANGINDOC) THEN C GET ANGLES FROM HEADER ANGBUF(1) = INUMBRT(K) CALL LUNGETVALS(INPROJ,IAPLOC+1,3,ANGBUF(2),IRTFLG) CALL BUILDM(INUMBRT,DM,1,ANGBUF,.FALSE.,SSDUM, & .FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ENDIF CLOSE(INPROJ) C PAD: PROJ TO SIZE: N CALL PADD2(PROJ,NS,BI,LSD,N) C FOURIER TRANSFORM OF: BI INV=+1 CALL FMRS_2(BI,N,N,INV) c$omp parallel do private(i,j) DO J=1,N DO I=0,N2 BI(I,J) = BI(I,J)*(-1)**(I+J+1) ENDDO ENDDO DO ISYM=1,MAXSYM IF (MAXSYM. GT. 1) THEN C SYMMETRIES, MULTIPLY MATRICES DMS = MATMUL(SM(:,:,ISYM),DM(:,:,K)) ELSE DMS = DM(:,:,K) ENDIF #ifdef SP_MP IF (RANDLIST(K)) THEN DO JT=1,LN1 C,schedule(static) c$omp parallel do private(j),shared(N,N2,JT,XE,WE,BI,DMS) DO J=-N2+JT,N2,LN1 CALL ONELINE(J,N,N2,XE,WE,BI,DMS) ENDDO ENDDO ELSE DO JT=1,LN1 c$omp parallel do private(j),shared(N,N2,JT,XO,WO,BI,DMS) DO J=-N2+JT,N2,LN1 CALL ONELINE(J,N,N2,XO,WO,BI,DMS) ENDDO ENDDO ENDIF #else IF (RANDLIST(K)) THEN DO J=-N2+1,N2 CALL ONELINE(J,N,N2,XE,WE,BI,DMS) ENDDO ELSE DO J=-N2+1,N2 CALL ONELINE(J,N,N2,XO,WO,BI,DMS) ENDDO ENDIF #endif C END OF SYMMETRIES LOOP ENDDO C END OF PROJECTIONS LOOP ENDDO #endif C SYMMETRIZE BOTH VOLUMES c$omp parallel sections c$omp section CALL SYMPLANE0(XE,WE,N2,N) c$omp section CALL SYMPLANE0(XO,WO,N2,N) c$omp end parallel sections 9999 IF (ALLOCATED(PROJ)) DEALLOCATE (PROJ) IF (ALLOCATED(BI)) DEALLOCATE (BI) IF (ALLOCATED(RANDLIST)) DEALLOCATE (RANDLIST) END C ------------------- WINDKB2 ------------------------------- SUBROUTINE WINDKB2(BI,R,L,LSD,N) DIMENSION R(L,L,L),BI(LSD,N,N) COMMON /BESSEL_PARAM/ ALPHA,AAAA,NNN PARAMETER (QUADPI = 3.14159265358979323846) PARAMETER (TWOPI = 2*QUADPI) IP = (N-L)/2+MOD(L,2) DO K=1,L DO J=1,L DO I=1,L R(I,J,K) = BI(IP+I,IP+J,IP+K) ENDDO ENDDO ENDDO L2 = (L/2)**2 L2P = (L/2-1)**2 IP = L/2+1 XNU = REAL(NNN)/2. RI = RIBSL(ALPHA,XNU) C IF (ABS(RI-RIN).GT.1.E-5) PRINT *,'BESSIK' WKB0 = ALPHA**XNU/RI QRT = (TWOPI*AAAA)**2 TNR = 0.0 M = 0 DO K=1,L DO J=1,L DO I=1,L LR = (K-IP)**2+(J-IP)**2+(I-IP)**2 IF (LR.LE.L2) THEN SIGMA=QRT*LR-ALPHA*ALPHA IF (ABS(SIGMA).LT.1.0E-7) THEN WKB=1.0 ELSEIF(SIGMA.GT.0.0) THEN C 2PI A R > ALPHA ART = SQRT(SIGMA) RI = RJBSL(ART, XNU) C if(abs(ri-rin)/rin.gt.1.e-5) print *,'bessjy',i,j,k WKB=WKB0*RI/ART**XNU ELSE C 2PI A R < ALPHA ART = SQRT(ABS(SIGMA)) RI = RIBSL(ART,XNU) C if(abs(ri-rin)/rin.gt.1.e-5) print *,'bessik',i,j,k,ri,rin WKB=WKB0*RI/ART**XNU ENDIF R(I,J,K) = R(I,J,K)/ABS(WKB) IF (LR.GE.L2P .AND. LR.LE.L2) THEN TNR=TNR+R(I,J,K) M=M+1 ENDIF ENDIF ENDDO ENDDO ENDDO TNR = TNR/REAL(M) c$omp parallel do private(i,j,k,lr) DO K=1,L DO J=1,L DO I=1,L LR=(K-IP)**2+(J-IP)**2+(I-IP)**2 IF (LR.LE.L2) THEN R(I,J,K)=R(I,J,K)-TNR ELSE R(I,J,K)=0.0 ENDIF ENDDO ENDDO ENDDO END C ----------------SYMPLANE0 --------------------------------------- SUBROUTINE SYMPLANE0(X,W,N2,N) DIMENSION W(0:N2,N,N) COMPLEX X(0:N2,N,N) C SYMMETRIZE PLANE 0 DO IZA=2,N2 DO IYA=2,N2 X(0,IYA,IZA)=X(0,IYA,IZA)+CONJG(X(0,N-IYA+2,N-IZA+2)) W(0,IYA,IZA)=W(0,IYA,IZA)+W(0,N-IYA+2,N-IZA+2) X(0,N-IYA+2,N-IZA+2)=CONJG(X(0,IYA,IZA)) W(0,N-IYA+2,N-IZA+2)=W(0,IYA,IZA) X(0,N-IYA+2,IZA)=X(0,N-IYA+2,IZA)+CONJG(X(0,IYA,N-IZA+2)) W(0,N-IYA+2,IZA)=W(0,N-IYA+2,IZA)+W(0,IYA,N-IZA+2) X(0,IYA,N-IZA+2)=CONJG(X(0,N-IYA+2,IZA)) W(0,IYA,N-IZA+2)=W(0,N-IYA+2,IZA) ENDDO ENDDO DO IYA=2,N2 X(0,IYA,1)=X(0,IYA,1)+CONJG(X(0,N-IYA+2,1)) W(0,IYA,1)=W(0,IYA,1)+W(0,N-IYA+2,1) X(0,N-IYA+2,1)=CONJG(X(0,IYA,1)) W(0,N-IYA+2,1)=W(0,IYA,1) ENDDO DO IZA=2,N2 X(0,1,IZA)=X(0,1,IZA)+CONJG(X(0,1,N-IZA+2)) W(0,1,IZA)=W(0,1,IZA)+W(0,1,N-IZA+2) X(0,1,N-IZA+2)=CONJG(X(0,1,IZA)) W(0,1,N-IZA+2)=W(0,1,IZA) ENDDO END C----------------ADDADA --------------------------------------- SUBROUTINE ADDADA(X,Y,N) DIMENSION X(N),Y(N) c$omp parallel do private(k) DO K=1,N X(K) = X(K)+Y(K) ENDDO END C++********************************************************************* C C WIW32D_DL.F 01/10/00 C ANGLES IN HEADER JULY 01 ArDean Leith C OPFILE 'U' FEB 02 ArDean Leith C BESSELS 05/03/02 P.A. Penczek C OPFILEC FEB 03 ArDean Leith C BUILDM PARAMETERS JUL 03 ArDean Leith C BUILDM PARAMETERS SEP 03 ArDean Leith C MPI OCT 03 CHAO YANG C REWRITE DEC 06 ArDean Leith C C=********************************************************************** C=* From: SPIDER - MODULAR IMAGE PROCESSING SYSTEM * C=* Copyright (C)2002, P. A. Penczek * C=* University of Texas - Houston Medical School * C=* Email: pawel.a.penczek@uth.tmc.edu * 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 WIW32D_DL(CALLRTSQ) C C23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 C--********************************************************************* SUBROUTINE WIW32D_DL(CALLRTSQ) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' INCLUDE 'F90ALLOC.INC' C DOC FILE POINTERS REAL, DIMENSION(:,:), POINTER :: ANGBUF, ANGSYM REAL, ALLOCATABLE, DIMENSION(:,:) :: DM,SM REAL, ALLOCATABLE, DIMENSION(:,:,:) :: WE,WO COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: XE,XO REAL, ALLOCATABLE, DIMENSION(:) :: TEMP REAL, ALLOCATABLE, DIMENSION(:) :: ROTANG,SHX,SHY LOGICAL :: ANGINDOC,CALLRTSQ CHARACTER(LEN=1) :: NULL CHARACTER(LEN=MAXNAM) :: FILPAT CHARACTER(LEN=MAXNAM) :: EVENVOL,ODDVOL,BOTHVOL CHARACTER(LEN=MAXNAM) :: ANGDOC C COMMON: /TABS/ IS USED IN ONELINE, EXTRACTLINE, PUTLINE3, ETC INTEGER, PARAMETER :: LTAB=4999 COMMON /TABS/ LN2,FLTB,TABI(0:LTAB) DATA IOPIC/98/,INPROJ/99/,LUNDOC/80/ #ifdef USE_MPI INCLUDE 'mpif.h' INTEGER MYPID, COMM, MPIERR COMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(COMM, MYPID, MPIERR) #else MYPID = -1 #endif NULL = CHAR(0) NILMAX = NIMAX CALL OPFILES(0,INPROJ,LUNDOC, FILPAT,NLET, 'O', & ITYPE,NSAM,NROW,NSLICE,MAXIM, & 'TEMPLATE FOR IMAGE FILES (E.G. STK@****)~', & .FALSE., INUMBR,NILMAX, NANG,IMGNUM, IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 MAXNUM = MAXVAL(INUMBR(1:NANG)) C NANG - TOTAL NUMBER OF IMAGES IF (MYPID .LE. 0) WRITE(NOUT,2001) NANG 2001 FORMAT(' NUMBER OF IMAGES: ',I7) C RETRIEVE ARRAY WITH ANGLES DATA IN IT ANGINDOC = .TRUE. C PSI, THE, PHI, REF#, EXP#, INPLANE, SX, SY MAXXT = 8 + 1 MAXYT = MAXNUM CALL GETDOCDAT('ANGLES DOC',.TRUE.,ANGDOC,77,.FALSE.,MAXXT, & MAXYT,ANGBUF,IRTFLG) IF (IRTFLG .NE. 0) ANGINDOC = .FALSE. C RETRIEVE ARRAY WITH SYMMETRIES DATA IN IT MAXXS = 0 MAXSYM = 0 CALL GETDOCDAT('SYMMETRIES DOC',.TRUE.,ANGDOC,77,.TRUE.,MAXXS, & MAXSYM,ANGSYM,IRTFLG) IF (IRTFLG .NE. 0) MAXSYM = 1 N2 = 2 * NSAM LSD = N2 + 2 - MOD(N2,2) NMAT = LSD * N2 * N2 NDIMANG = 1 IF (ANGINDOC) NDIMANG = MAXNUM ALLOCATE(DM(9,NDIMANG), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'BP 32F, DM',9 * NDIMANG) GOTO 999 ENDIF IF (ANGINDOC) THEN C GET ANGLES FROM DOCUMENT FILE CALL BUILDM1(INUMBR,DM,9,NANG,ANGBUF,.FALSE.,SSDUM, & .TRUE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 C NO LONGER NEED ANGBUF IF (.NOT. CALLRTSQ) DEALLOCATE(ANGBUF) ENDIF NDIMSYM = MAX(1,MAXSYM) ALLOCATE(SM(9,NDIMSYM), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'BP 32F, SM',9*NDIMSYM) GOTO 999 ENDIF IF (MAXSYM .GT. 1) THEN C HAVE SYMMETRIES CALL BUILDS(SM,MAXSYM,ANGSYM(1,1),IRTFLG) DEALLOCATE(ANGSYM) ENDIF ALLOCATE(XE(0:NSAM, N2,N2), & WE(0:NSAM, N2,N2), & XO(0:NSAM, N2,N2), & WO(0:NSAM, N2,N2), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN C SOME ARRAYS ARE COMPLEX SO 6 NOT 4 MEMWANT = 6 * (NSAM+1) * N2 * N2 CALL ERRT(46,'BP 32F; XE, WE, WO, & XO',MEMWANT) GOTO 999 ENDIF C CREATE FFTW3 PLAN FOR 3D FFT ON XE USING ALL THREADS CALL FMRS_PLAN(.TRUE.,XE,N2,N2,N2, 0,-1,IRTFLG) IF (IRTFLG .NE. 0) RETURN C GENERALIZED KAISER-BESSEL WINDOW ACCORDING TO LEWITT LN = 5 LN2 = LN / 2 ! USED IN COMMON /TABS/ R = NSAM / 2 V = REAL(LN-1) / 2.0 / REAL(N2) ALPHA = 6.5 AAAA = 0.9*V NNN = 3 C GENERATE TABLE WITH INTERPOLANTS B0 = SQRT(ALPHA) * BESI1(ALPHA) FLTB = REAL(LTAB) / REAL(LN2+1) cc$omp parallel do private(i,s,x),shared(mmm) DO I=0,LTAB S = REAL(I) / FLTB / N2 IF (S .LE. AAAA) THEN X = SQRT(1.0 - (S/AAAA)**2) TABI(I) = SQRT(ALPHA*X) * BESI1(ALPHA*X) / B0 ELSE TABI(I) = 0.0 ENDIF ENDDO C WIW32DQ(NS,XE,WE,XO,WO,LSD,N,N2, C WIW32DQ_DL(NSAM,XE,WE,XO,WO,LSD,N CALL WIW32DQ_DL(NSAM, XE,WE,XO,WO,LSD,N2, CALLRTSQ,FILPAT, & INPROJ,ANGBUF,INUMBR,NANG, & DM,IMGNUM,SM,MAXSYM,ANGINDOC,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL FILERD(BOTHVOL,NLETI,NULL,'RECONSTRUCTED VOLUME',IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL FILERD(EVENVOL,NLETI,NULL,'FIRST SAMPLE VOLUME',IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL FILERD(ODDVOL,NLETI,NULL,'SECOND SAMPLE VOLUME',IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 C CBE AND CBO (XE & XO) ARE LSD x N2 x N2 C CWE AND CWO (WE & WO) ARE (LSD/2) x N2 x N2 LSDD2 = LSD / 2 C STORE FOURIER XE TEMPORARILY IN OVERALL VOLUME: BOTHVOL MAXIM = 0 ITYPE = 3 CALL OPFILEC(0,.FALSE.,BOTHVOL,IOPIC,'U',ITYPE,LSD,N2,N2, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL WRITEV(IOPIC,XE,LSD,N2,LSD,N2,N2) IF (MYPID .LE. 0) CLOSE(IOPIC) C STORE WEIGHT WE TEMPORARILY IN ODDVOL (ODD VOLUME) MAXIM = 0 CALL OPFILEC(0,.FALSE.,ODDVOL,IOPIC,'U',ITYPE,LSDD2,N2,N2, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL WRITEV(IOPIC,WE,LSDD2,N2,LSDD2,N2,N2) IF (MYPID .LE. 0) CLOSE(IOPIC) C PROCESS EVEN VOLUME -------------------------------- EVEN c write(6,*) 'vals: even',xe(1,1,1),xe(1,1,1), c & we(1,1,1),we(2,2,2), c & xo(2,2,2),xo(2,2,2), c & wo(2,2,2),wo(2,2,2) C WEIGHT, FOURIER TRANSFORM, & WINDOW USING: XE AND WE CALL NRMW2(XE,WE,NSAM,N2) CALL WINDKB2A(XE,XE,NSAM,LSD,N2,ALPHA,AAAA,NNN) C NOW EVEN VOLUME IS READY, SYMMETRIZE IF NECESSARY IF (MAXSYM .GT. 1) THEN C ADDITIONAL SYMMETRIZATION OF VOLUME XE IN REAL SPACE 05/03/02 ALLOCATE (TEMP(NSAM*NSAM*NSAM), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN MEMWANT = NSAM * NSAM * NSAM CALL ERRT(46,'BP 32F, TEMP',MEMWANT) GOTO 999 ENDIF CALL COP(XE,TEMP,NSAM*NSAM*NSAM) c$omp parallel do private(i,j,k) DO K=1,N2 DO J=1,N2 DO I=0,NSAM XE(I,J,K) = CMPLX(0.0,0.0) ENDDO ENDDO ENDDO IF (MOD(NSAM,2) .EQ. 0) THEN KNX = NSAM/2-1 ELSE KNX = NSAM/2 ENDIF KLX = -NSAM/2 CALL SYMVOL(TEMP,XE,KLX,KNX,KLX,KNX,KLX,KNX,SM,MAXSYM) ENDIF C WRITE EVEN VOLUME TO FILE: EVENVOL, NOTE: NSAM=NROW=NSLICE MAXIM = 0 CALL OPFILEC(0,.FALSE.,EVENVOL,IOPIC,'U',ITYPE,NSAM,NSAM,NSAM, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL WRTVOL(IOPIC,NSAM,NSAM, 1,NSAM, XE,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 IF (MYPID .LE. 0) CLOSE(IOPIC) C RECOVER XE FROM BOTHVOL AND PUT IT IN: XE --------------- BOTH MAXIM = 0 CALL OPFILEC(0,.FALSE.,BOTHVOL,IOPIC,'O',ITYPE,LSD,N2,N2, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL READV(IOPIC,XE,LSD,N2,LSD,N2,N2) CLOSE(IOPIC) C RECOVER WEIGHT WE TEMP. STORED IN: ODDVOL AND PUT IT IN: WE MAXIM = 0 CALL OPFILEC(0,.FALSE.,ODDVOL,IOPIC,'O',ITYPE,LSDD2,N2,N2, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL READV(IOPIC,WE,LSDD2,N2,LSDD2,N2,N2) CLOSE(IOPIC) c write(6,*) 'vals: rec ',xe(1,1,1),xe(1,1,1), c & we(1,1,1),we(2,2,2), c & xo(2,2,2),xo(2,2,2), c & wo(2,2,2),wo(2,2,2) C ADD EVEN AND ODD VOLUMES, PLACE IN XE AND WE CALL ADDADA(XE,XO,NMAT) CALL ADDADA(WE,WO,NMAT/2) C WEIGHT, FOURIER TRANSFORM, AND WINDOW BOTHVOL USING: XE AND WE CALL NRMW2(XE,WE,NSAM,N2) CALL WINDKB2A(XE,XE,NSAM,LSD,N2,ALPHA,AAAA,NNN) C OVEALL VOLUME IS NOW IN: XE IF (MAXSYM .GT. 1) THEN C ADDITIONAL SYMMETRIZATION OF VOLUME TOTAL IN REAL SPACE 05/03/02 CALL COP(XE, TEMP,NSAM*NSAM*NSAM) c$omp parallel do private(i,j,k) DO K=1,N2 DO J=1,N2 DO I=0,NSAM XE(I,J,K) = CMPLX(0.0,0.0) ENDDO ENDDO ENDDO CALL SYMVOL(TEMP,XE,KLX,KNX,KLX,KNX,KLX,KNX,SM,MAXSYM) ENDIF C WRITE OVERALL VOLUME FROM: XE TO: BOTHVOL MAXIM = 0 CALL OPFILEC(0,.FALSE.,BOTHVOL,IOPIC,'U',ITYPE,NSAM,NSAM,NSAM, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL WRTVOL(IOPIC,NSAM,NSAM, 1,NSAM, XE,IRTFLG) IF (MYPID .LE. 0) CLOSE(IOPIC) c write(6,*) 'vals: end ',xe(1,1,1),xe(1,1,1), c & we(1,1,1),we(2,2,2), c & xo(2,2,2),xo(2,2,2), c & wo(2,2,2),wo(2,2,2) C ODD VOLUME IS STILL UNDAMAGED IN: XO --------------------- ODD C WEIGHT, FOURIER TRANSFORM, AND WINDOW ODD VOL. USING: XO AND WO CALL NRMW2( XO,WO,NSAM,N2) CALL WINDKB2A(XO,XO,NSAM,LSD,N2,ALPHA,AAAA,NNN) IF (MAXSYM .GT. 1) THEN C ADDITIONAL SYMMETRIZATION OF VOLUME XO IN REAL SPACE 05/03/02 CALL COP(XO,TEMP, NSAM*NSAM*NSAM) c$omp parallel do private(i,j,k) DO K=1,N2 DO J=1,N2 DO I=0,NSAM XO(I,J,K) = CMPLX(0.0,0.0) ENDDO ENDDO ENDDO CALL SYMVOL(TEMP,XO,KLX,KNX,KLX,KNX,KLX,KNX,SM,MAXSYM) ENDIF C WRITE ODD FILE VOLUME FROM: XO TO: ODDVOL MAXIM = 0 CALL OPFILEC(0,.FALSE.,ODDVOL,IOPIC,'U',ITYPE,NSAM,NSAM,NSAM, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL WRTVOL(IOPIC,NSAM,NSAM, 1,NSAM, XO,IRTFLG) IF (MYPID .LE. 0) CLOSE(IOPIC) 999 CONTINUE IF (ALLOCATED(DM)) DEALLOCATE(DM) IF (ALLOCATED(SM)) DEALLOCATE(SM) IF (ALLOCATED(XE)) DEALLOCATE(XE) IF (ALLOCATED(WE)) DEALLOCATE(WE) IF (ALLOCATED(XO)) DEALLOCATE(XO) IF (ALLOCATED(WO)) DEALLOCATE(WO) IF (ALLOCATED(TEMP)) DEALLOCATE(TEMP) IF (ASSOCIATED(ANGBUF)) DEALLOCATE(ANGBUF) IF (ASSOCIATED(ANGSYM)) DEALLOCATE(ANGSYM) CLOSE(INPROJ) END C ---------------- WIW32DQ_DL ----------------------------------- SUBROUTINE WIW32DQ_DL(NSAM,XE,WE,XO,WO,LSD,N, CALLRTSQ,FILPAT, & INPROJ,ANGBUF,INUMBRT,NANG, & DM,IMGNUM,SM,MAXSYM,ANGINDOC,IRTFLG) C NOTE: STUPID TRANSFORM OF N2-->N !!!!al INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' REAL :: DM(3,3,NANG),SM(3,3,MAXSYM),DMS(3,3) REAL, DIMENSION(9,NANG) :: ANGBUF LOGICAL :: ANGINDOC,CALLRTSQ LOGICAL, ALLOCATABLE, DIMENSION(:) :: RANDLIST REAL, ALLOCATABLE, DIMENSION(:,:) :: PROJ,PROJTEMP COMPLEX, ALLOCATABLE, DIMENSION(:,:) :: BI REAL, DIMENSION(4) :: ANGBUFT INTEGER, DIMENSION(NANG) :: INUMBRT REAL :: WE(0:NSAM,N,N) REAL :: WO(0:NSAM,N,N) COMPLEX :: XE(0:NSAM,N,N) COMPLEX :: XO(0:NSAM,N,N) CHARACTER(LEN=MAXNAM) :: FILNAM,ODDVOL,FILPAT CHARACTER (LEN=MAXNAM) :: FILPATOUT INTEGER, ALLOCATABLE, DIMENSION(:) :: INUMBROUT DOUBLE PRECISION :: PI LOGICAL :: ITMP DATA LUNROTT/81/,LUNDOC/80/ #ifdef USE_MPI c This MPI version is memory intensive. c It requires 6 copies of the 3-D volume, c 2-D images are read into memory and distributed. c Each processor will hold roughly nang/nproc 2-D images. c Memory requirement may be reduced in a future version? INCLUDE 'mpif.h' INTEGER MYPID, COMM, MPIERR, NPROCS, K3 INTEGER ISTAT(MPI_STATUS_SIZE) INTEGER, ALLOCATABLE, DIMENSION(:) :: PSIZE, NBASE INTEGER NANGLOC, NREM, IPROC, ISAM, JROW, JGLB, JLOC REAL , ALLOCATABLE, DIMENSION(:,:,:) :: PRJLOC REAL , ALLOCATABLE, DIMENSION(:,:,:) :: PRJBUF REAL , ALLOCATABLE, DIMENSION(:,:,:) :: WELOC, WOLOC COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: XELOC, XOLOC COMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(COMM, MYPID, MPIERR) CALL MPI_COMM_SIZE(COMM, NPROCS, MPIERR) #else MYPID = -1 #endif c$omp parallel do private(i,j,k) DO K=1,N DO J=1,N DO I=0,NSAM XE(I,J,K) = CMPLX(0.0,0.0) WE(I,J,K) = 0.0 XO(I,J,K) = CMPLX(0.0,0.0) WO(I,J,K) = 0.0 ENDDO ENDDO ENDDO NSIZE = 1 LUNROT = 0 ! NO RTSQ OUTPUT WANTED IF (CALLRTSQ) THEN NSIZE = NSAM ALLOCATE(INUMBROUT(NANG), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'INUMBROUT', NANG) GOTO 9999 ENDIF C OPEN OUTPUT IMAGE(S) FOR RTSQ ITYPE = 1 CALL OPFILES(INPROJ,LUNROTT,LUNDOC, FILPATOUT,NLETO, 'U', & ITYPE,NSAM,NSAM,1, MAXIM, & 'TRANSFORMED OUTPUT IMAGES TEMPLATE (E.G. ROT@****)~', & .FALSE., INUMBROUT,NANG, NANGT, IMGNUMOUT, IRTFLG) IF (IRTFLG .EQ. 0) LUNROT = LUNROTT IF (IRTFLG .GT. 0) GOTO 9999 ENDIF ALLOCATE (PROJ(NSAM,NSAM), & PROJTEMP(NSAM,NSIZE), & BI(0:NSAM,N), & RANDLIST(NANG), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN MEMWANT = NSAM*NSAM + NSAM*NSIZE + 2*(NSAM+1)*N + NANG CALL ERRT(46,'BP 32F, PROJ, ...',MEMWANT) RETURN ENDIF C CREATE FFTW3 PLAN FOR 2D FFT ON BI USING ALL THREADS CALL FMRS_PLAN(.TRUE.,BI,N,N,1, 0,+1,IRTFLG) IF (IRTFLG .NE. 0) RETURN C CREATE LIST OF IMAGES FOR EACH RECONSTRUCTION RANDLIST(1:NANG/2) = .TRUE. RANDLIST(NANG/2+1:NANG) = .FALSE. DO K=1,NANG CALL RANDOM_NUMBER(HARVEST=X) IORD = MIN(NANG,MAX(1,INT(X*NANG+0.5))) ITMP = RANDLIST(IORD) RANDLIST(IORD) = RANDLIST(K) RANDLIST(K) = ITMP ENDDO #ifdef USE_MPI C --- BEGIN MPI VERSION ----------------------------------- MPI C DISTRIBUTE PARTICLES TO PROCESSORS. C NANGLOC IS THE NUMBER OF PARTICLES ASSIGNED TO EACH PROCESSOR. ALLOCATE(PSIZE(NPROCS) ,NBASE(NPROCS), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MEMWANT = 2 * NPROCS CALL ERRT(46,'BP 32F, PSIZE & NBASE',MEMWANT) RETURN ENDIF C SETPART RETURNS THE SIZE OF THE LOCAL PIECE AND GLOBAL OFFSET. CALL SETPART(NANG, PSIZE, NBASE) NANGLOC = PSIZE(MYPID+1) C 2-D IMAGES ARE DISTRIBUTED AND HELD IN PRJLOC ON EACH PROCESSOR ALLOCATE(PRJBUF(NSAM,NSAM,PSIZE(1)), & PRJLOC(NSAM,NSAM,NANGLOC), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MEMWANT = NSAM*NSAM*PSIZE(1) + NSAM*NSAM*NANGLOC CALL ERRT(46,'BP 32F, PRJBUF, PRJLOC',MEMWANT) RETURN ENDIF C PROCESSOR 0 READS IMAGE FILES AND DISTRIBUTE THEM. C (THIS VERSION ASSUMES THAT THERE IS SUFFICIENT C MEMORY TO HOLD NANG/NPROCS IMAGES) DO IPROC = 1, NPROCS NANGLOC = PSIZE(IPROC) C READ IMAGES INTO THE BUFFER FIRST, THEN DISTRIBUTED DO JLOC = 1, NANGLOC JGLB = NBASE(IPROC) + JLOC CALL FILGET(FILPAT,FILNAM,NLET,INUMBR(JGLB),IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 MAXIM = 0 CALL OPFILEC(0 ,.FALSE., FILNAM, INPROJ, 'O' , & IFORM , NSAM , NSAM , NSL , MAXIM, & 'DUMMY',.FALSE., IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 CALL READV1P(INPROJ,PRJBUF(1,1,JLOC), & NSAM,NSAM,NSAM,NSAM,1) CLOSE(INPROJ) ENDDO IF (IPROC .GT. 1) THEN IF (MYPID .EQ. 0) THEN C SEND TO ANOTHER PROCESSOR #ifdef MPI_DEBUG WRITE(6,*) 'WIW32DQ: SENDING TO PID = ', IPROC-1 #endif CALL MPI_SEND(PRJBUF , NSAM*NSAM*NANGLOC, MPI_REAL, & IPROC-1, IPROC-1 , COMM , & MPIERR) IF (MPIERR .NE. 0) THEN WRITE(6,*) 'WIW32DQ: SEND ERROR!' STOP ENDIF ELSEIF (MYPID .EQ. IPROC-1) THEN C RECEIVE PROJECTION IMAGES FROM PROCESSOR 0 CALL MPI_RECV(PRJLOC, NSAM*NSAM*NANGLOC, MPI_REAL, & 0 , MPI_ANY_TAG , COMM , & ISTAT , MPIERR) IF (MPIERR .NE. 0) THEN WRITE(6,*) 'WIW32DQ: RECV FAILED' STOP ENDIF #ifdef MPI_DEBUG WRITE(6,111) MYPID 111 FORMAT(1X,'WIW32DQ: MYPID = ', I3, ' RECEIVED!') #endif ENDIF ELSEIF (MYPID .EQ. 0) THEN C KEEP FOR MYSELF DO JLOC = 1, NANGLOC DO ISAM = 1, NSAM DO JROW = 1, NSAM PRJLOC(ISAM,JROW,JLOC) = PRJBUF(ISAM,JROW,JLOC) ENDDO ENDDO ENDDO ENDIF ENDDO IF (ALLOCATED(PRJBUF)) DEALLOCATE(PRJBUF) C IF (.NOT. ANGINDOC) THEN C GET ANGLES FROM HEADER ANGBUFT(1) = INUMBRT(K) CALL LUNGETVALS(INPROJ,IAPLOC+1,8,ANGBUFT(2),IRTFLG) CALL BUILDM(INUMBRT,DM,1,ANGBUFT,.FALSE.,SSDUM, & .FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ENDIF CLOSE(INPROJ) C PERFORM CALCULATIONS IN PARALLEL NOW ALLOCATE (WELOC(0:NSAM,N,N), XELOC(0:NSAM,N,N), & WOLOC(0:NSAM,N,N), XOLOC(0:NSAM,N,N), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MEMWANT = 4*(NSAM+1)*N*N CALL ERRT(46,'BP 3F, WELOC, XELOC, ...',MEMWANT) RETURN ENDIF DO K=1,N DO J=1,N DO I=0,NSAM XELOC(I,J,K) = CMPLX(0.0,0.0) XOLOC(I,J,K) = CMPLX(0.0,0.0) WELOC(I,J,K) = 0.0 WOLOC(I,J,K) = 0.0 ENDDO ENDDO ENDDO C NANGLOC = PSIZE(MYPID+1) DO JLOC = 1, NANGLOC JGLB = NBASE(MYPID+1) + JLOC C PAD: PRJLOC TO SIZE: N CALL PADD2(PRJLOC(1,1,JLOC),NSAM,BI,LSD,N) C FOURIER TRANSFORM OF: BI INV = +1 CALL FMRS_2(BI,N,N,INV) DO J=1,N DO I=0,NSAM BI(I,J)=BI(I,J)*(-1)**(I+J+1) ENDDO ENDDO DO ISYM=1,MAXSYM IF (MAXSYM .GT. 1 ) THEN C SYMMETRIES, MULTIPLY MATRICES DMS = MATMUL(SM(:,:,ISYM),DM(:,:,JGLB)) ELSE DMS = DM(:,:,JGLB) ENDIF IF (RANDLIST(JGLB)) THEN DO J=-NSAM+1,NSAM CALL ONELINE(J,N,NSAM,XELOC,WELOC,BI,DMS) ENDDO ELSE DO J=-NSAM+1,NSAM CALL ONELINE(J,N,NSAM,XOLOC,WOLOC,BI,DMS) ENDDO ENDIF ENDDO ENDDO IF (ALLOCATED(PRJLOC)) DEALLOCATE(PRJLOC) C SUM UP X AND W FROM THE LOCAL PIECES (XLOC, WLOC) C RESIDING ON EACH PROCESSOR DO K3 = 1, N CALL MPI_ALLREDUCE(XELOC(1,1,K3), XE(1,1,K3), (NSAM+1)*N, & MPI_COMPLEX , MPI_SUM , COMM , & MPIERR) IF ( MPIERR. NE. 0 ) THEN WRITE(6,*) 'WIW32DQ: FAILED TO ALLREDUCE' STOP ENDIF CALL MPI_ALLREDUCE(XOLOC(1,1,K3), XO(1,1,K3), (NSAM+1)*N, & MPI_COMPLEX , MPI_SUM , COMM , & MPIERR) IF ( MPIERR. NE. 0 ) THEN WRITE(6,*) 'WIW32DQ: FAILED TO ALLREDUCE' STOP ENDIF CALL MPI_ALLREDUCE(WELOC(1,1,K3), WE(1,1,K3), (NSAM+1)*N, & MPI_REAL , MPI_SUM , COMM , & MPIERR) IF ( MPIERR. NE. 0 ) THEN WRITE(6,*) 'WIW32DQ: FAILED TO ALLREDUCE' STOP ENDIF CALL MPI_ALLREDUCE(WOLOC(1,1,K3), WO(1,1,K3), (NSAM+1)*N, & MPI_REAL , MPI_SUM , COMM , & MPIERR) IF ( MPIERR. NE. 0 ) THEN WRITE(6,*) 'WIW32DQ: FAILED TO ALLREDUCE' STOP ENDIF ENDDO IF (ALLOCATED(XOLOC)) DEALLOCATE (XOLOC) IF (ALLOCATED(XELOC)) DEALLOCATE (XELOC) IF (ALLOCATED(WOLOC)) DEALLOCATE (WOLOC) IF (ALLOCATED(WELOC)) DEALLOCATE (WELOC) #ifdef MPI_DEBUG WRITE(6,*) 'WIW32DQ: COMPLETED GLOBAL SUM, MYPID = ', MYPID #endif C -------------------- END OF MPI VERSION ---------------------- #else #ifdef SP_MP LN1 = LN + 1 ! WHY?? LN ALWAYS EQUALED: 5 (IN CALLER) LN1 = 6 #endif NWANT = 1 DO IF (CALLRTSQ) THEN C READ IMAGE AND ROTATE, SCALE & SHIFT INTO: PROJ CALL REDVOL(INPROJ,NSAM,NSAM, 1,1,PROJTEMP,IRTFLG) C Reg. numbers for angle, scale,& shift =(6,0,7,8) IF (VERBOSE) WRITE(NOUT,90)IMGNUM, & ANGBUF(7,IMGNUM), & ANGBUF(8,IMGNUM), & ANGBUF(9,IMGNUM) 90 FORMAT(' IMAGE: ',I6, & ' ANGLE: ',G10.3, & ' SHIFT: (',G10.3,',',G10.3,')') CALL ROT2QS_DL(PROJTEMP,PROJ, NSAM,NSAM, & ANGBUF(7,IMGNUM),1.0,ANGBUF(8,IMGNUM), & ANGBUF(9,IMGNUM),0,0) IF (LUNROT .GT. 0) THEN CALL WRTVOL(LUNROT,NSAM,NSAM,1,1, PROJ,IRTFLG) ENDIF ELSE C READ IMAGE INTO: PROJ CALL REDVOL(INPROJ,NSAM,NSAM, 1,1, PROJ,IRTFLG) ENDIF IF (.NOT. ANGINDOC) THEN C GET PROJECTION ANGLES FROM HEADER ANGBUFT(1) = INUMBRT(NWANT) CALL LUNGETVALS(INPROJ,IAPLOC+1,8,ANGBUFT(2),IRTFLG) CALL BUILDM(INUMBRT,DM,1,ANGBUFT,.FALSE.,SSDUM, & .FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ENDIF C PAD: PROJ TO SIZE: N CALL PADD2(PROJ,NSAM,BI,LSD,N) C FOURIER TRANSFORM OF: BI INV = +1 CALL FMRS_2(BI,N,N,INV) c$omp parallel do private(i,j) DO J=1,N DO I=0,NSAM BI(I,J) = BI(I,J) * (-1)**(I+J+1) ENDDO ENDDO DO ISYM=1,MAXSYM IF (MAXSYM. GT. 1) THEN C SYMMETRIES, MULTIPLY MATRICES DMS = MATMUL(SM(:,:,ISYM),DM(:,:,IMGNUM)) ELSE DMS = DM(:,:,IMGNUM) ENDIF #ifdef SP_MP IF (RANDLIST(NWANT)) THEN DO JT=1,LN1 c$omp parallel do private(j) DO J=-NSAM+JT,NSAM,LN1 CALL ONELINE(J,N,NSAM,XE,WE,BI,DMS) ENDDO ENDDO ELSE DO JT=1,LN1 c$omp parallel do private(j) DO J=-NSAM+JT,NSAM,LN1 CALL ONELINE(J,N,NSAM,XO,WO,BI,DMS) ENDDO ENDDO ENDIF #else IF (RANDLIST(NWANT)) THEN DO J=-NSAM+1,NSAM CALL ONELINE(J,N,NSAM,XE,WE,BI,DMS) ENDDO ELSE DO J=-NSAM+1,NSAM CALL ONELINE(J,N,NSAM,XO,WO,BI,DMS) ENDDO ENDIF #endif ENDDO ! END OF SYMMETRIES LOOP IF (NWANT .GE. NANG) EXIT ! END OF INPUT LIST C OPEN NEXT SET OF I/O FILES CALL NEXTFILES(NWANT, INUMBR,INUMBROUT, & INPROJ,INPROJ,LUNROT, & FILPAT,FILPATOUT, & IMGNUM,IMGNUMOUT,IRTFLG) IF (IRTFLG .LT. 0) EXIT ! END OF INPUT STACK IF (IRTFLG .NE. 0) GOTO 9999 ! ERROR ENDDO ! END OF PROJECTIONS LOOP #endif C SYMMETRIZE PLANE 0 OF BOTH VOLUMES c$omp parallel sections c$omp section CALL SYMPLANE0(XE,WE,NSAM,N) c$omp section CALL SYMPLANE0(XO,WO,NSAM,N) c$omp end parallel sections 9999 IF (ALLOCATED(PROJ)) DEALLOCATE(PROJ) IF (ALLOCATED(BI)) DEALLOCATE(BI) IF (ALLOCATED(RANDLIST)) DEALLOCATE(RANDLIST) IF (ALLOCATED(INUMBROUT))DEALLOCATE(INUMBROUT) CLOSE(LUNROTT) END C ------------------- WINDKB2A ------------------------------- SUBROUTINE WINDKB2A(BI,R,L,LSD,N,ALPHA,AAAA,NNN) DIMENSION R(L,L,L), BI(LSD,N,N) PARAMETER (QUADPI = 3.14159265358979323846) PARAMETER (TWOPI = 2*QUADPI) IP = (N-L )/ 2 + MOD(L,2) DO K=1,L DO J=1,L DO I=1,L R(I,J,K) = BI(IP+I,IP+J,IP+K) ENDDO ENDDO ENDDO L2 = (L/2)**2 L2P = (L/2-1)**2 IP = L / 2+1 XNU = REAL(NNN) / 2. RI = RIBSL(ALPHA,XNU) WKB0 = ALPHA**XNU / RI QRT = (TWOPI*AAAA)**2 TNR = 0.0 M = 0 DO K=1,L DO J=1,L DO I=1,L LR = (K-IP)**2+(J-IP)**2+(I-IP)**2 IF (LR.LE.L2) THEN SIGMA = QRT*LR-ALPHA*ALPHA IF (ABS(SIGMA).LT.1.0E-7) THEN WKB=1.0 ELSEIF(SIGMA.GT.0.0) THEN C 2PI A R > ALPHA ART = SQRT(SIGMA) RI = RJBSL(ART, XNU) WKB = WKB0*RI/ART**XNU ELSE C 2PI A R < ALPHA ART = SQRT(ABS(SIGMA)) RI = RIBSL(ART,XNU) WKB = WKB0*RI/ART**XNU ENDIF R(I,J,K) = R(I,J,K) / ABS(WKB) IF (LR.GE.L2P .AND. LR.LE.L2) THEN TNR = TNR+R(I,J,K) M = M+1 ENDIF ENDIF ENDDO ENDDO ENDDO TNR = TNR / REAL(M) c$omp parallel do private(i,j,k,lr) DO K=1,L DO J=1,L DO I=1,L LR = (K-IP)**2+(J-IP)**2+(I-IP)**2 IF (LR .LE. L2) THEN R(I,J,K) = R(I,J,K) - TNR ELSE R(I,J,K) = 0.0 ENDIF ENDDO ENDDO ENDDO END