C++********************************************************************* C C WIW3D.F 05/03/02 C REPLACED BESSEL FUNCTIONS C OPFILEC FEB 03 ARDEAN LEITH C BUILDM PARAMETERS JUL 03 ARDEAN LEITH C BUILDM PARAMETERS SEP 03 ARDEAN LEITH C ALLOCATION ERROR HANDLING OCT 04 ARDEAN LEITH C OMP PRIVATE XX NOV 04 ARDEAN LEITH C MPI DECONSTRUCTED FROM WIW32D NOV 06 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 WIW3D C C PURPOSE: FOR 'BP 3F' BACK PROJECTION - 3D, INTERPOLATED IN C FOURIER SPACE ||* C AND FOR: 'BPD 3F' IMPROVED WITH LESS MEMORY REQUIRED AND C MORE STACK USAGE. C C23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 C--********************************************************************* SUBROUTINE WIW3D INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' INCLUDE 'F90ALLOC.INC' REAL, DIMENSION(:,:), ALLOCATABLE :: DM,SM COMPLEX, DIMENSION(:,:,:), ALLOCATABLE :: X REAL, DIMENSION(:), ALLOCATABLE :: TEMP C DOC FILE POINTERS REAL, DIMENSION(:,:), POINTER :: ANGBUF, ANGSYM COMMON /F_SPEC/ FINPAT,NLET,FINPIC CHARACTER*80 FINPIC,FINPAT,FILNAM,ANGDOC DATA IOPIC/98/,INPIC/99/ #ifdef USE_MPI C KLUDGE TO OVERCOME BUG REPORTED BY L. ALAMO, NOV 2006 CALL WIW3D_MPI() RETURN #else MYPID = -1 #endif 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 N - LINEAR DIMENSION OF PROJECTIONS AND RESTORED CUBE C NANG - TOTAL NUMBER OF IMAGES IF (MYPID .LE. 0) WRITE(NOUT,2001) NANG 2001 FORMAT(' NUMBER OF IMAGES =',I5) C RETRIEVE ARRAY WITH ANGLES DATA IN IT MAXXT = 4 MAXYT = MAXNUM CALL GETDOCDAT('ANGLES DOC',.TRUE.,ANGDOC,77,.FALSE.,MAXXT, & MAXYT,ANGBUF,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 C RETRIEVE ARRAY WITH SYMMETRIES DATA IN IT MAXXS = 4 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),INTFLG) MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINPIC,INPIC,'O',IFORM,NSAM,NROW,NSL, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN CLOSE(INPIC) N2 = 2*NSAM LSD = N2+2-MOD(N2,2) NMAT = LSD*N2*N2 ALLOCATE(DM(9,NANG), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'BP 3F, DM',9*NANG) GOTO 999 ENDIF CALL BUILDM(INUMBR,DM,NANG,ANGBUF(1,1),.FALSE.,SSDUM, & .FALSE.,IRTFLG) DEALLOCATE(ANGBUF) IF (IRTFLG .NE. 0) GOTO 999 IF (MAXSYM .GT. 1) THEN C HAS SYMMETRIES ALLOCATE(SM(9,MAXSYM), & X(0:N2/2,N2,N2), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,' BP 3F, SM...',9*MAXSYM+NSAM*N2*N2) GOTO 999 ENDIF CALL BUILDS(SM,MAXSYM,ANGSYM(1,1),IRTFLG) DEALLOCATE(ANGSYM) ELSE ALLOCATE(SM(1,1), & X(0:N2/2,N2,N2), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'BP 3F, SM...',1+NSAM*N2*N2) GOTO 999 ENDIF ENDIF C CREATE FFTW3 PLAN FOR 3D FFT ON X USING ALL THREADS CALL FMRS_PLAN(.TRUE.,X,N2,N2,N2, 0,-1,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL WIW3DQ(NSAM,X,LSD,N2,N2/2,INUMBR,DM,NANG,SM,MAXSYM,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 C ADDITIONAL SYMMETRIZATION OF THE VOLUME IN REAL SPACE 05/03/02 IF (MAXSYM .GT. 1) THEN ALLOCATE(TEMP(NSAM*NSAM*NSAM), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'BP 3F, TEMP',NSAM*NSAM*NSAM) GOTO 999 ENDIF CALL COP(X,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 X(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,X,KLX,KNX,KLX,KNX,KLX,KNX,SM,MAXSYM) DEALLOCATE(TEMP) ENDIF IFORM = 3 CALL OPFILEC(0,.TRUE.,FILNAM,IOPIC,'U',IFORM,NSAM,NSAM,NSAM, & MAXIM,'RECONSTRUCTED 3-D',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9997 C NOTE: NSAM=NROW=NSLICE CALL WRITEV(IOPIC,X,NSAM,NSAM,NSAM,NSAM,NSAM) 9997 CLOSE(IOPIC) 999 IF (ALLOCATED(DM)) DEALLOCATE(DM) IF (ALLOCATED(SM)) DEALLOCATE(SM) IF (ALLOCATED(X)) DEALLOCATE(X) END C ------------------ WIW3DQ ---------------------------------- SUBROUTINE WIW3DQ(NS,X,LSD,N,N2,INUMBR, & DM,NANG,SM,MAXSYM,IRTFLG) C NOTE: STUPID TRANSFORM OF N2-->N AND N2/2-->N2 !!!!al COMPLEX X(0:N2,N,N) DIMENSION INUMBR(NANG) DIMENSION DM(3,3,NANG),SM(3,3,MAXSYM),DMS(3,3) REAL, DIMENSION(:,:), ALLOCATABLE :: PROJ COMPLEX, DIMENSION(:,:), ALLOCATABLE :: BI REAL, DIMENSION(:,:,:), ALLOCATABLE :: W CHARACTER*80 FINPIC,FINPAT COMMON /F_SPEC/ FINPAT,NLET,FINPIC DOUBLE PRECISION PI PARAMETER (LTAB=4999) COMMON /TABS/ LN2,FLTB,TABI(0:LTAB) COMMON /BESSEL_PARAM/ ALPHA,AAAA,NNN C,mmm DATA INPROJ/99/ PARAMETER (QUADPI = 3.1415926535897932384626) PARAMETER (TWOPI = 2*QUADPI) #ifdef USE_MPI c This MPI version is memory intensive. c It requires 4 copies of the 3-D volumn, c 2-D images are read into memory and distributed. c Each processor will hold roughly nang/nproc c 2-D images. Memory requirement will be reduced c in a future version. INCLUDE 'mpif.h' INTEGER MYPID, COMM, MPIERR, NPROCS INTEGER ISTAT(MPI_STATUS_SIZE) INTEGER, ALLOCATABLE, DIMENSION(:) :: PSIZE, NBASE INTEGER NANGLOC, NREM, IPROC, ISAM, JROW, JGLB, JLOC REAL , ALLOCATABLE, DIMENSION(:,:,:) :: PRJLOC, PRJBUF REAL , ALLOCATABLE, DIMENSION(:,:,:) :: WLOC COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: XLOC COMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(COMM, MYPID , MPIERR) CALL MPI_COMM_SIZE(COMM, NPROCS, MPIERR) #else MYPID = -1 #endif LN = 5 LN2 = LN/2 C Generalized Kaiser-Bessel window according to Lewitt R = NS/2 V = REAL(LN-1)/2.0/REAL(N) 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,xx) DO I=0,LTAB S=REAL(I)/FLTB/N IF (S.LE.AAAA) THEN XX = SQRT(1.0-(S/AAAA)**2) TABI(I) = SQRT(ALPHA*XX)*BESI1(ALPHA*XX)/B0 ELSE TABI(I) = 0.0 ENDIF ENDDO ALLOCATE(W(0:N2,N,N), & PROJ(NS,NS), & BI(0:N2,N), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN NEEDED = (N2+1)*N*N + NS*NS + 2*(N2+1)*N CALL ERRT(46,'BP 3F, W,PROJ,BI',NEEDED) RETURN ENDIF C CREATE FFTW3 PLAN FOR 2D FFT ON BI USING ALL THREADS CALL FMRS_PLAN(.TRUE.,BI,N2,N2,1, 0,+1,IRTFLG) IF (IRTFLG .NE. 0) RETURN c$omp parallel do private(i,j,k) DO K=1,N DO J=1,N DO I=0,N2 X(I,J,K) = CMPLX(0.0,0.0) W(I,J,K) = 0.0 ENDDO ENDDO ENDDO #ifdef SP_MP LN1=LN+1 #endif #ifdef USE_MPI C DISTRIBUTE PARTICLES TO PROCESSORS. C NANGLOC IS THE NUMBER OF PARTICLES ASSIGNED C TO EACH PROCESSOR. ALLOCATE(PSIZE(NPROCS),NBASE(NPROCS),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'BP 3F, PSIZE,NBASE',NPROCS*2) RETURN ENDIF C SETPART RETURNS THE SIZE OF THE LOCAL PIECE C AND THE GLOBAL OFFSET. CALL SETPART(NANG, PSIZE, NBASE) NANGLOC = PSIZE(MYPID+1) C 2-D IMAGES ARE DISTRIBUTED AND HELD IN PRJLOC C ON EACH PROCESSOR ALLOCATE(PRJBUF(NS,NS,PSIZE(1)), PRJLOC(NS,NS,NANGLOC), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN NEEDED = NS*NS*PSIZE(1) + NS*NS*NANGLOC CALL ERRT(46,'BP 3F, PRJLOC...',NEEDED) 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 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 999 MAXIM = 0 CALL OPFILEC(0 ,.FALSE., FINPIC, INPROJ, 'O' , & IFORM , NSAM , NSAM , NSL , MAXIM, & 'DUMMY',.FALSE., IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 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,*) 'WIW3DQ: 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,*) 'WIW3DQ: 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,*) 'WIW3DQ: RECV FAILED' STOP END IF #ifdef MPI_DEBUG WRITE(6,111) MYPID 111 FORMAT(' WIW3DQ: 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) END DO END DO END DO ENDIF ENDDO IF (ALLOCATED(PRJBUF)) DEALLOCATE(PRJBUF) C PERFORM CALCULATIONS IN PARALLEL NOW ALLOCATE (WLOC(0:N2,N,N), XLOC(0:N2,N,N), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = ((N2+1)*N*N) * 2 CALL ERRT(46,'BP 3F, WLOC, XLOC',MWANT) RETURN ENDIF DO K=1,N DO J=1,N DO I=0,N2 XLOC(I,J,K)=CMPLX(0.0,0.0) WLOC(I,J,K)=0.0 ENDDO ENDDO ENDDO NANGLOC = PSIZE(MYPID+1) DO JLOC = 1, NANGLOC JGLB = NBASE(MYPID+1) + JLOC C PAD: BI 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 DO J=-N2+1,N2 CALL ONELINE(J,N,N2,XLOC,WLOC,BI,DMS) ENDDO END DO ENDDO IF (ALLOCATED(PRJLOC)) DEALLOCATE(PRJLOC) C SUM UP X AND W FROM THE LOCAL PIECES (XLOC, WLOC) C RESIDING ON EACH PROCESSOR DO K = 1, N CALL MPI_ALLREDUCE(XLOC(1,1,K), X , (N2+1)*N, MPI_COMPLEX, & MPI_SUM , COMM, MPIERR) IF ( MPIERR. NE. 0 ) THEN WRITE(6,*) 'WIW3DQ: FAILED TO ALLREDUCE' STOP END IF CALL MPI_ALLREDUCE(WLOC(1,1,K), W , (N2+1)*N, MPI_REAL, & MPI_SUM , COMM, MPIERR) IF ( MPIERR. NE. 0 ) THEN WRITE(6,*) 'WIW3DQ: FAILED TO ALLREDUCE' STOP END IF ENDDO IF (ALLOCATED(XLOC)) DEALLOCATE (XLOC) IF (ALLOCATED(WLOC)) DEALLOCATE (WLOC) #ifdef MPI_DEBUG WRITE(6,*) 'WIW3DQ: COMPLETED GLOBAL SUM, MYPID = ', MYPID #endif #else DO K=1,NANG C PRINT *,' PROJECTION #',K C OPEN DESIRED FILE CALL FILGET(FINPAT,FINPIC,NLET,INUMBR(K),IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINPIC,INPROJ,'O',IFORM,NSAM,NSAM,NSL, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL READV(INPROJ,PROJ,NS,NS,NS,NS,1) CLOSE(INPROJ) C PAD: PROJ TO SIZE: N CALL PADD2(PROJ,NS,BI,LSD,N) C FOURIER TRANSFOR 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 DO JT=1,LN1 c$omp parallel do private(j) DO J=-N2+JT,N2,LN1 CALL ONELINE(J,N,N2,X,W,BI,DMS) ENDDO ENDDO #else DO J=-N2+1,N2 CALL ONELINE(J,N,N2,X,W,BI,DMS) ENDDO #endif ENDDO ! END OF SYMMETRIES LOOP ENDDO ! END OF PROJECTIONS LOOP #endif C SYMMETRIZE PLANE 0 CALL SYMPLANE0(X,W,N2,N) CALL NRMW2(X,W,N2,N) CALL WINDKB2(X,X,NS,LSD,N) IRTFLG = 0 999 CONTINUE IF (ALLOCATED(W)) DEALLOCATE (W) IF (ALLOCATED(PROJ)) DEALLOCATE (PROJ) IF (ALLOCATED(BI)) DEALLOCATE (BI) END C ------------------- PADD2 ------------------------------- SUBROUTINE PADD2(PROJ,L,BI,LSD,N) C PADS: PROJ OF SIZE: L INTO: BI WITH SIZE: N DIMENSION PROJ(L,L),BI(LSD,N) DOUBLE PRECISION QS KLP = 0 R = L/2 QS = 0.0D0 C ESTIMATE AVERAGE OUTSIDE THE CIRCLE CALL ASTA(PROJ,L,R,QS,KLP) QS = QS/REAL(KLP) C ZEROS ALL OF: BI c$omp parallel do private(i,j) DO J=1,N DO I=1,N BI(I,J) = 0.0 ENDDO ENDDO C FOR L ODD ADD ONE. N IS ALWAYS EVEN IP = (N-L)/2+MOD(L,2) c$omp parallel do private(i,j) DO J=1,L DO I=1,L BI(IP+I,IP+J) = PROJ(I,J) - QS ENDDO ENDDO END C ------------------- NRMW2 ------------------------------- SUBROUTINE NRMW2(R,W,N2,N) DIMENSION W(0:N2,N,N) COMPLEX R(0:N2,N,N) c$omp parallel do private(i,j,k) DO K=1,N DO J=1,N DO I=0,N2 IF (W(I,J,K).GT.0.1) THEN R(I,J,K) = R(I,J,K)*(-1)**(I+J+K)/W(I,J,K) ELSE R(I,J,K) = (0.0,0.0) ENDIF ENDDO ENDDO ENDDO INV = -1 CALL FMRS_3(R,N,N,N,INV) END C++********************************************************************* C C WIW3D_MPI.F DECONSTRUCTED FROM WIW32D NOV 06 ARDEAN LEITH C ORIGINAL WIW3D GAVE WRONG RESULTS C WHEN USED UNDER MPI C C=********************************************************************** C=* From: SPIDER - MODULAR IMAGE PROCESSING SYSTEM * C=* Copyright (C)2002, 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 WIW3D_MPI C C23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 C--********************************************************************* SUBROUTINE WIW3D_MPI 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 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: WE REAL, DIMENSION(:), ALLOCATABLE :: TEMP LOGICAL :: ANGINDOC CHARACTER(LEN=1) :: NULL CHARACTER(LEN=MAXNAM) :: ANGDOC CHARACTER(LEN=MAXNAM) :: FINPAT,FINPIC CHARACTER(LEN=MAXNAM) :: FOUTPIC #ifndef SP_32 INTEGER *8 :: FFTW_PLAN_F,FFTW_PLAN_B #else INTEGER *4 :: FFTW_PLAN_F,FFTW_PLAN_B #endif DATA IOPIC/18/,INPIC/19/ #ifdef USE_MPI INCLUDE 'mpif.h' INTEGER MYPID, COMM, MPIERR COMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(COMM, MYPID, MPIERR) #else PRINT *, ' THIS ROUTINE for MPI ONLY' STOP #endif #ifdef USE_MPI 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,0,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), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN C X ARRAYS ARE COMPLEX SO 3 NOT 2 MEMWANT = 3 * ((N2/2)+1) * N2 * N2 CALL ERRT(46,'BP 32F; XE & WE ',MEMWANT) GOTO 999 ENDIF C CREATE FFTW3 PLAN FOR 3D FFT ON X USING ALL THREADS CALL FMRS_PLAN(.TRUE.,XE,N2,N2,N2, 0,-1,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL WIW3DQ_MPI(NSAM,XE,WE,LSD,N2,N2/2,FINPAT,FINPIC,INPIC, & INUMBR,DM,NANG,SM,MAXSYM,ANGINDOC,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 CALL FILERD(FOUTPAC,NLETI,NULL,'RECONSTRUCTED 3D',IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 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 C OUTPUT FILE IFORM = 3 MAXIM = 0 CALL OPFILEC(0,.FALSE.,FOUTPIC,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) 999 CONTINUE IF (ALLOCATED(DM)) DEALLOCATE(DM) IF (ALLOCATED(SM)) DEALLOCATE(SM) IF (ALLOCATED(XE)) DEALLOCATE(XE) IF (ALLOCATED(WE)) DEALLOCATE(WE) IF (ALLOCATED(TEMP)) DEALLOCATE(TEMP) IF (ASSOCIATED(ANGBUF)) DEALLOCATE(ANGBUF) IF (ASSOCIATED(ANGSYM)) DEALLOCATE(ANGSYM) #endif END C ---------------- WIW3DQ_MPI ------------------------------------- SUBROUTINE WIW3DQ_MPI(NS,XE,WE,LSD,N,N2,FINPAT,FINPIC,INPROJ, & INUMBRT,DM,NANG,SM,MAXSYM,ANGINDOC,IRTFLG) C NOTE: STUPID TRANSFORM OF N2-->N AND N2/2-->N2 !!!!al INCLUDE 'CMLIMIT.INC' LOGICAL :: ANGINDOC REAL, DIMENSION(:,:), ALLOCATABLE :: PROJ COMPLEX, DIMENSION(:,:), ALLOCATABLE :: BI REAL, DIMENSION(4) :: ANGBUF DIMENSION :: WE(0:N2,N,N) COMPLEX :: XE(0:N2,N,N) INTEGER :: INUMBRT(NANG) REAL :: DM(3,3,NANG) REAL :: SM(3,3,MAXSYM),DMS(3,3) CHARACTER(LEN=*) :: FINPAT,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 PARAMETER (QUADPI = 3.1415926535897932384) PARAMETER (TWOPI = 2*QUADPI) c This MPI version is memory intensive. c 2-D images are read into memory and distributed. c Each processor will hold roughly nang/nproc c 2-D images. #ifndef USE_MPI MYPID = -1 PRINT *, ' THIS ROUTINE FOR MPI COMPILATION/USE ONLY' STOP #else INCLUDE 'mpif.h' INTEGER MYPID, COMM, MPIERR, NPROCS, K3 INTEGER ISTAT(MPI_STATUS_SIZE) INTEGER, ALLOCATABLE, DIMENSION(:) :: PSIZE, NBASE REAL , ALLOCATABLE, DIMENSION(:,:,:) :: PRJLOC REAL , ALLOCATABLE, DIMENSION(:,:,:) :: PRJBUF REAL , ALLOCATABLE, DIMENSION(:,:,:) :: WELOC COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: XELOC COMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(COMM, MYPID, MPIERR) CALL MPI_COMM_SIZE(COMM, NPROCS, MPIERR) LN = 5 LN2 = LN/2 C GENERALIZED KAISER-BESSEL WINDOW ACCORDING TO LEWITT R = NS/2 V = REAL(LN-1)/2.0/REAL(N) 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) 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 ENDDO ENDDO ENDDO ALLOCATE (PROJ(NS,NS), & BI(0:N2,N), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN MEMWANT = NS*NS + 2*(N2 + 1) * N 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,N2,N2,1, 0,+1,IRTFLG) IF (IRTFLG .NE. 0) RETURN 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 C AND THE 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 DISTRIBUTE DO JLOC = 1, NANGLOC JGLB = NBASE(IPROC) + JLOC CALL FILGET(FINPAT,FINPIC,0,INUMBRT(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 ELSEIF (MYPID .EQ. IPROC-1) THEN C C RECEIVE PROJECTION IMAGES FROM PROCESSOR 0 C 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 ELSEIF (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), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MEMWANT = 2*(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) WELOC(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 DO J=-N2+1,N2 CALL ONELINE(J,N,N2,XELOC,WELOC,BI,DMS) ENDDO 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(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 ENDDO IF (ALLOCATED(XELOC)) DEALLOCATE (XELOC) IF (ALLOCATED(WELOC)) DEALLOCATE (WELOC) #ifdef MPI_DEBUG WRITE(6,*) 'WIW32DQ: COMPLETED GLOBAL SUM, MYPID = ', MYPID #endif C --- END OF MPI VERSION --- C SYMMETRIZE VOLUME CALL SYMPLANE0(XE,WE,N2,N) 9999 IF (ALLOCATED(PROJ)) DEALLOCATE (PROJ) IF (ALLOCATED(BI)) DEALLOCATE (BI) #endif END C++********************************************************************* C C WIW3D_DL.F 05/03/02 C REPLACED BESSEL FUNCTIONS C OPFILEC FEB 03 ARDEAN LEITH C BUILDM PARAMETERS JUL 03 ARDEAN LEITH C BUILDM PARAMETERS SEP 03 ARDEAN LEITH C ALLOCATION ERROR HANDLING OCT 04 ARDEAN LEITH C OMP PRIVATE XX NOV 04 ARDEAN LEITH C MPI DECONSTRUCTED FROM WIW32D NOV 06 ARDEAN LEITH C REWRITE DEC 06 ArDean Leith 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 WIW3D_DL(CALLRTSQ) C C23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 C--********************************************************************* SUBROUTINE WIW3D_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 COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: X REAL, ALLOCATABLE, DIMENSION(:) :: TEMP REAL, ALLOCATABLE, DIMENSION(:) :: ROTANG,SHX,SHY LOGICAL :: ANGINDOC,CALLRTSQ CHARACTER(LEN=MAXNAM) :: FILPAT,VOLNAM CHARACTER(LEN=MAXNAM) :: ANGDOC DATA IOPIC/98/,INPIC/99/,LUNDOC/80/ #ifdef USE_MPI C KLUDGE TO OVERCOME BUG REPORTED BY L. ALAMO, NOV 2006 CALL WIW3D_MPI_DL() RETURN #else MYPID = -1 #endif NILMAX = NIMAX CALL OPFILES(0,INPIC,LUNDOC, FILPAT,NLET, 'O', & ITYPE,NSAM,NROW,NSLICE,MAXIM, & 'IMAGE FILE NAME OR TEMPLATE (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,I7) C RETRIEVE ARRAY WITH ANGLES DATA IN IT 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) GOTO 999 C RETRIEVE ARRAY WITH SYMMETRIES DATA IN IT MAXXS = 4 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 ALLOCATE(DM(9,MAXNUM), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'BP 3F, DM', 9*NANG) GOTO 999 ENDIF C GET ANGLES FROM DOCUMENT FILE AND PLACE IN DM CALL BUILDM1(INUMBR,DM,9,NANG,ANGBUF,.FALSE.,SSDUM, & .TRUE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 IF (.NOT. CALLRTSQ) DEALLOCATE(ANGBUF) NDIMSYM = MAX(1,MAXSYM) ALLOCATE(SM(9,NDIMSYM), X(0:NSAM,N2,N2), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN MWANT = 9*NDIMSYM + (NSAM+1)*N2*N2 CALL ERRT(46,'BP 3F, SM & X',MWANT) GOTO 999 ENDIF IF (MAXSYM .GT. 1) THEN C HAVE SYMMETRIES, CONSTRUCT SM ANGLES ARRAY CALL BUILDS(SM,MAXSYM,ANGSYM(1,1),IRTFLG) DEALLOCATE(ANGSYM) ENDIF CALL WIW3DQ_DL(NSAM,X,LSD,N2, CALLRTSQ,FILPAT,ANGBUF, & INUMBR,NANG, DM,IMGNUM,SM,MAXSYM,IRTFLG) IF (IRTFLG .NE. 0) GOTO 999 C write(6,*) ' 0 x(65,65,65):', x(65,65,65) C ADDITIONAL SYMMETRIZATION OF THE VOLUME IN REAL SPACE 05/03/02 IF (MAXSYM .GT. 1) THEN ALLOCATE(TEMP(NSAM*NSAM*NSAM), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'BP 3F, TEMP',NSAM*NSAM*NSAM) GOTO 999 ENDIF CALL COP(X,TEMP, NSAM*NSAM*NSAM) c$omp parallel do private(i,j,k) DO K=1,N2 DO J=1,N2 DO I=0,NSAM X(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,X,KLX,KNX,KLX,KNX,KLX,KNX,SM,MAXSYM) DEALLOCATE(TEMP) ENDIF IFORM = 3 CALL OPFILEC(0,.TRUE.,VOLNAM,IOPIC,'U',IFORM,NSAM,NSAM,NSAM, & MAXIM,'RECONSTRUCTED 3-D',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9997 C NOTE: NSAM=NROW=NSLICE CALL WRTVOL(IOPIC,NSAM,NSAM, 1,NSAM, X,IRTFLG) 9997 CLOSE(IOPIC) 999 IF (ALLOCATED(DM)) DEALLOCATE(DM) IF (ALLOCATED(SM)) DEALLOCATE(SM) IF (ALLOCATED(X)) DEALLOCATE(X) END C ------------------ WIW3DQ ---------------------------------- SUBROUTINE WIW3DQ_DL(NSAM,X,LSD,N2, CALLRTSQ,FILPAT,ANGBUF, & INUMBRT,NANG, DM,IMGNUM,SM,MAXSYM,IRTFLG) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' REAL :: DM(3,3,NANG),SM(3,3,MAXSYM),DMS(3,3) COMPLEX :: X(0:NSAM, N2,N2) CHARACTER(LEN=*) :: FILPAT REAL, DIMENSION(9,NANG) :: ANGBUF INTEGER, DIMENSION(NANG) :: INUMBRT REAL, ALLOCATABLE, DIMENSION(:,:) :: PROJ,PROJTEMP REAL, ALLOCATABLE, DIMENSION(:,:,:) :: W COMPLEX, ALLOCATABLE, DIMENSION(:,:) :: BI CHARACTER (LEN=MAXNAM) :: FILPATOUT,DOCNAM INTEGER, ALLOCATABLE, DIMENSION(:) :: INUMBROUT LOGICAL :: ANGINDOC,CALLRTSQ DOUBLE PRECISION :: PI C COMMON: /TABS/ IS USED IN ONELINE, EXTRACTLINE, PUTLINE3, ETC INTEGER, PARAMETER :: LTAB=4999 COMMON /TABS/ LN2,FLTB,TABI(0:LTAB) DATA INPROJ/99/,LUNROTT/81/,LUNDOC/80/ MYPID = -1 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 XT = SQRT(1.0 - (S/AAAA)**2) TABI(I) = SQRT(ALPHA*XT) * BESI1(ALPHA*XT) / B0 ELSE TABI(I) = 0.0 ENDIF 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 (IF ANY) 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) C IF IRTFLG IS NEGATIVE THEN NO RTSQ OUTPUT WANTED IF (IRTFLG .EQ. 0) LUNROT = LUNROTT IF (IRTFLG .GT. 0) GOTO 9999 ENDIF ALLOCATE(PROJ(NSAM,NSAM), & PROJTEMP(NSAM,NSIZE), & W(0:NSAM,N2,N2), & BI(0:NSAM,N2), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN NEEDED = NSAM*NSAM + NSAM*NSIZE + 2*(NSAM+1)*N2*N2 + & (NSAM+1)*N2 CALL ERRT(46,'BP 3F, PROJ...',NEEDED) RETURN ENDIF C CREATE FFTW3 PLAN FOR 2D FFT ON BI USING ALL THREADS CALL FMRS_PLAN(.TRUE.,BI,N2,N2,1, 0,+1,IRTFLG) IF (IRTFLG .NE. 0) RETURN c$omp parallel do private(i,j,k) DO K=1,N2 DO J=1,N2 DO I=0,NSAM X(I,J,K) = CMPLX(0.0,0.0) W(I,J,K) = 0.0 ENDDO ENDDO ENDDO #ifdef SP_MP C LN1 = LN + 1 ! WHY?? LN ALWAYS EQUALS: 5 LN1 = 6 ! WHY?? LN ALWAYS EQUALS: 5 #endif NWANT = 1 DO c write(6,*) 'Projection #: ',IMGNUM IF (CALLRTSQ) THEN C READ IMAGE AND ROTATE, SCALE & SHIFT INTO: PROJ CALL REDVOL(INPROJ,NSAM,NSAM, 1,1,PROJTEMP,IRTFLG) 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,')') C Reg. numbers for angle, scale,& shift =(6,0,7,8) 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 c write(6,*) ' 0 maxval(proj):', maxval(proj) C PAD: PROJ TO SIZE: N2 CALL PADD2(PROJ,NSAM,BI,LSD,N2) C FOURIER TRANSFORM OF: BI INV = +1 CALL FMRS_2(BI,N2,N2,INV) c$omp parallel do private(i,j) DO J=1,N2 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 DO JT=1,LN1 c$omp parallel do private(j) DO J=-NSAM+JT,NSAM,LN1 CALL ONELINE(J,N2,NSAM,X,W,BI,DMS) ENDDO ENDDO #else DO J=-NSAM+1,NSAM CALL ONELINE(J,N2,NSAM,X,W,BI,DMS) ENDDO #endif ENDDO ! END OF SYMMETRIES LOOP IF (NWANT .GE. NANG) EXIT ! END OF LIST C OPEN NEXT SET OF I/O FILES CALL NEXTFILES(NWANT, INUMBRT,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 c write(6,*) ' 0 x(65,65,65):', x(65,65,65) C SYMMETRIZE PLANE 0 CALL SYMPLANE0(X,W,NSAM,N2) C WEIGHT, FOURIER TRANSFORM, AND WINDOW CALL NRMW2(X,W,NSAM,N2) CALL WINDKB2A(X,X,NSAM,LSD,N2,ALPHA,AAAA,NNN) IRTFLG = 0 9999 IF (ALLOCATED(W)) DEALLOCATE (W) IF (ALLOCATED(PROJ)) DEALLOCATE (PROJ) IF (ALLOCATED(BI)) DEALLOCATE (BI) CLOSE(LUNROTT) END C++********************************************************************* C C WIW3D_MPI.F DECONSTRUCTED FROM WIW32D NOV 06 ARDEAN LEITH C ORIGINAL WIW3D GAVE WRONG RESULTS C WHEN USED UNDER MPI C C ********************************************************************** C C WIW3D_MPI C C23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 C--********************************************************************* SUBROUTINE WIW3D_MPI_DL C NOTE: STUPID TRANSFORM OF N2-->N AND N2/2-->N2 !!!!al INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' INCLUDE 'F90ALLOC.INC' C DOC FILE POINTERS REAL, DIMENSION(:,:), POINTER :: ANGBUF, ANGSYM REAL, ALLOCATABLE, DIMENSION(:,:) :: DM,SM COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: XE REAL, ALLOCATABLE, DIMENSION(:,:,:) :: WE REAL, ALLOCATABLE, DIMENSION(:) :: TEMP LOGICAL :: ANGINDOC,CALLRTSQ CHARACTER(LEN=1) :: NULL CHARACTER(LEN=MAXNAM) :: ANGDOC CHARACTER(LEN=MAXNAM) :: FILPAT,FILNAM CHARACTER(LEN=MAXNAM) :: FOUTPIC DATA IOPIC/18/,INPIC/19/ #ifdef USE_MPI INCLUDE 'mpif.h' INTEGER MYPID, COMM, MPIERR COMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(COMM, MYPID, MPIERR) #else PRINT *, ' THIS ROUTINE for MPI ONLY' STOP #endif #ifdef USE_MPI NULL = CHAR(0) NILMAX = NIMAX CALL FILELIST(.TRUE.,INPIC,FILPAT,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. 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 C OPEN FIRST IMAGE FILE TO DETERMINE NSAM, NROW, NSL CALL FILGET(FILPAT,FILNAM,0,INUMBR(1),IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 MAXIM = 0 CALL OPFILEC(0,.FALSE.,FILNAM,INPIC,'O',IFORM,NSAM,NROW,NSL, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 CLOSE(INPIC) N2 = 2 * NSAM LSD = N2 + 2 - MOD(N2,2) NMAT = LSD * N2 * N2 IF (ANGINDOC) THEN C GET ANGLES FROM DOC. FILE ALLOCATE(DM(9,MAXNUM), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN MEMWANT = 9 * NANG CALL ERRT(46,'BP 32F, DM',MEMWANT) GOTO 9999 ENDIF CALL BUILDM1(INUMBR,DM,9,NANG,ANGBUF,.FALSE.,SSDUM, & .TRUE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 IF (ASSOCIATED(ANGBUF)) DEALLOCATE(ANGBUF) ELSE C GET ANGLES FROM IMAGE FILE HEADER ALLOCATE(DM(9,1), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'BP 32F, DM',IER) GOTO 9999 ENDIF ENDIF IF (MAXSYM .GT. 1) THEN C HAVE SYMMETRIES ALLOCATE(SM(9,MAXSYM), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN MEMWANT = 9 * MAXSYM CALL ERRT(46,'BP 3F, SM',MEMWANT) GOTO 9999 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 9999 ENDIF ENDIF ALLOCATE(XE(0:NSAM,N2,N2), & WE(0:NSAM,N2,N2), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN C X ARRAYS ARE COMPLEX SO 3 NOT 2 MEMWANT = 3 * ((NSAM)+1) * N2 * N2 CALL ERRT(46,'BP 32F; XE & WE ',MEMWANT) GOTO 9999 ENDIF CALL WIW3DQ_MPI_DL(NSAM,XE,WE,LSD,N2, FILPAT,FILNAM,INPIC, & INUMBR,NANG, DM,NANG,SM,MAXSYM,ANGINDOC,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 c write(6,*) ' x4(65,65,65):', x(65,65,65) CALL FILERD(FOUTPAC,NLETI,NULL,'RECONSTRUCTED 3D',IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C WEIGHT, FOURIER TRANSFORM, AND WINDOW CALL NRMW2(XE,WE,NSAM,N2) CALL WINDKB2A(XE,XE,NSAM,LSD,N2,ALPHA,AAAA,NNN) 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 9999 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 OPEN OUTPUT VOLUME IFORM = 3 MAXIM = 0 CALL OPFILEC(0,.FALSE.,FOUTPIC,IOPIC,'U',IFORM,NSAM,NSAM,NSAM, & MAXIM,' ',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C NOTE: NSAM=NROW=NSLICE CALL WRTVOL(IOPIC,NSAM,NSAM, 1,NSAM, XE,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 IF (MYPID .LE. 0) CLOSE(IOPIC) 9999 IF (ALLOCATED(DM)) DEALLOCATE(DM) IF (ALLOCATED(SM)) DEALLOCATE(SM) IF (ALLOCATED(XE)) DEALLOCATE(XE) IF (ALLOCATED(WE)) DEALLOCATE(WE) IF (ALLOCATED(TEMP)) DEALLOCATE(TEMP) IF (ASSOCIATED(ANGBUF)) DEALLOCATE(ANGBUF) IF (ASSOCIATED(ANGSYM)) DEALLOCATE(ANGSYM) #endif END C ---------------- WIW3DQ_MPI ------------------------------------- SUBROUTINE WIW3DQ_MPI_DL(NSAM,XE,WE,LSD,N2, & FILPAT,FILNAM,INPROJ, & INUMBRT,NGOT,DM, NANG,SM,MAXSYM,ANGINDOC,IRTFLG) INCLUDE 'CMLIMIT.INC' LOGICAL :: ANGINDOC REAL, DIMENSION(:,:), ALLOCATABLE :: PROJ COMPLEX, DIMENSION(:,:), ALLOCATABLE :: BI REAL, DIMENSION(4) :: ANGBUF DIMENSION :: WE(0:NSAM,N2,N2) COMPLEX :: XE(0:NSAM,N2,N2) INTEGER :: INUMBRT(NANG) REAL :: DM(3,3,NANG) REAL :: SM(3,3,MAXSYM),DMS(3,3) CHARACTER(LEN=*) :: FILPAT,FILNAM DOUBLE PRECISION :: PI LOGICAL :: ITMP C COMMON: /TABS/ IS USED IN ONELINE, EXTRACTLINE, PUTLINE3, ETC PARAMETER (LTAB=4999) COMMON /TABS/ LN2,FLTB,TABI(0:LTAB) c This MPI version is memory intensive. c 2-D images are read into memory and distributed. c Each processor will hold roughly nang/nproc c 2-D images. #ifndef USE_MPI MYPID = -1 PRINT *, ' THIS ROUTINE FOR MPI COMPILATION AND USE ONLY' STOP #else INCLUDE 'mpif.h' INTEGER MYPID, COMM, MPIERR, NPROCS, K3 INTEGER ISTAT(MPI_STATUS_SIZE) INTEGER, ALLOCATABLE, DIMENSION(:) :: PSIZE, NBASE REAL , ALLOCATABLE, DIMENSION(:,:,:) :: PRJLOC REAL , ALLOCATABLE, DIMENSION(:,:,:) :: PRJBUF REAL , ALLOCATABLE, DIMENSION(:,:,:) :: WELOC COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: XELOC COMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(COMM, MYPID, MPIERR) CALL MPI_COMM_SIZE(COMM, NPROCS, MPIERR) LN = 5 LN2 = LN/2 C GENERALIZED KAISER-BESSEL WINDOW ACCORDING TO LEWITT 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) 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 DO K=1,N2 DO J=1,N2 DO I=0,NSAM XE(I,J,K) = CMPLX(0.0,0.0) WE(I,J,K) = 0.0 ENDDO ENDDO ENDDO ALLOCATE (PROJ(NSAM,NSAM), & BI(0:NSAM,N2), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN MEMWANT = NSAM*NSAM + (NSAM + 1) * N2 CALL ERRT(46,'BP 32F, PROJ, ...',MEMWANT) RETURN ENDIF 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 C AND THE 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 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 DISTRIBUTE DO JLOC = 1, NANGLOC JGLB = NBASE(IPROC) + JLOC CALL FILGET(FILPAT,FILNAM,0,INUMBRT(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,*)'WIW32DQ: MYPID=', MYPID,' 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) IF (.NOT. ANGINDOC) THEN C GET ANGLES FROM HEADER ANGBUF(1) = INUMBRT(K) CALL LUNGETVALS(INPROJ,IAPLOC+1,3,ANGBUF(2),IRTFLG) CALL BUILDM1(INUMBRT,DM,4,1,ANGBUF,.FALSE.,SSDUM, & .TRUE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ENDIF CLOSE(INPROJ) C PERFORM CALCULATIONS IN PARALLEL NOW ALLOCATE (WELOC(0:NSAM,N2,N2), XELOC(0:NSAM,N2,N2), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MEMWANT = 2 * (NSAM+1)*N2*N2 CALL ERRT(46,'BP 3F, WELOC, XELOC',MEMWANT) RETURN ENDIF DO K=1,N2 DO J=1,N2 DO I=0,NSAM XELOC(I,J,K) = CMPLX(0.0,0.0) WELOC(I,J,K) = 0.0 ENDDO ENDDO ENDDO NANGLOC = PSIZE(MYPID+1) DO JLOC = 1, NANGLOC JGLB = NBASE(MYPID+1) + JLOC C PAD: PRJLOC TO SIZE: N2 CALL PADD2(PRJLOC(1,1,JLOC),NSAM,BI,LSD,N2) C FOURIER TRANSORM OF: BI INV = +1 CALL FMRS_2(BI,N2,N2,INV) DO J=1,N2 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 DO J=-NSAM+1,NSAM CALL ONELINE(J,N2,NSAM,XELOC,WELOC,BI,DMS) ENDDO 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)*N2, & 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)*N2, & MPI_REAL , MPI_SUM , COMM , & MPIERR) IF ( MPIERR. NE. 0 ) THEN WRITE(6,*) 'WIW32DQ: FAILED TO ALLREDUCE' STOP ENDIF ENDDO IF (ALLOCATED(XELOC)) DEALLOCATE (XELOC) IF (ALLOCATED(WELOC)) DEALLOCATE (WELOC) #ifdef MPI_DEBUG WRITE(6,*) 'WIW32DQ: COMPLETED GLOBAL SUM, MYPID = ', MYPID #endif C SYMMETRIZE VOLUME CALL SYMPLANE0(XE,WE,NSAM,N2) 9999 IF (ALLOCATED(PROJ)) DEALLOCATE(PROJ) IF (ALLOCATED(BI)) DEALLOCATE(BI) #endif END C++********************************************************************* C C BUILDM1.F MERGED WITH REANG JUL 03 ARDEAN LEITH C BYLIST ADDED SEP 03 ARDEAN LEITH C ********************************************************************** C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM. AUTHOR: J.FRANK * C=* Copyright (C) 1985-2005 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 BUILDM1(ILIST,DM,IX,NANG,ANGBUF,FILLSS,SS,BYLIST,IRTFLG) C C PURPOSE: BULID ROTATION MATRICES FROM THREE EULERIAN ANGLES C C23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 C--********************************************************************* SUBROUTINE BUILDM1(ILIST,DM,IX,NANG,ANGBUF,FILLSS, & SS,BYLIST,IRTFLG) INCLUDE 'CMBLOCK.INC' REAL :: DM(9,NANG),ANGBUF(IX,NANG),SS(6,NANG) INTEGER :: ILIST(NANG) LOGICAL :: FILLSS,BYLIST #ifdef USE_MPI include 'mpif.h' icomm = MPI_COMM_WORLD call MPI_COMM_RANK(icomm, mypid, ierr) #else mypid = -1 #endif C READ ANGLES FROM THE DOCUMENT FILE. C ORDER IN THE DOCUMENT FILE IS PSI, THETA, PHI AND ANGLES C ARE IN DEGREES! IN ANG ARRAY IT IS OTHER WAY AROUND C OUTPUT IS COMPACTED TO 1...NANG LINES (NOT BY SELECTOR) DO K=1,NANG C GET ANGLE SELECTOR INDEX FROM ILIST ITMP = ILIST(K) ICOUNT = ANGBUF(1,ITMP) IF (ICOUNT .LE. 0) THEN C MISSING KEY CALL ERRT(102,'MISSING ANGLE FOR IMAGE',ITMP) IRTFLG = 1 RETURN ENDIF KT = K IF (BYLIST) KT = ITMP CALL CANG(ANGBUF(4,ITMP),ANGBUF(3,ITMP),ANGBUF(2,ITMP), & FILLSS,SS(1,KT),DM(1,KT)) IF (VERBOSE) THEN IF (MYPID .LE. 0)WRITE(NOUT,333)K,(ANGBUF(J,ITMP),J=2,4) 333 FORMAT(' PROJECTION #',I7, & '; PSI=',F6.1,' THETA=',F6.1,' PHI=',F6.1) ENDIF ENDDO IRTFLG = 0 END