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 MPI CHANGES OCT 06 ArDean Leith C MPI READV1P BUG DEC 09 ArDean Leith C C=********************************************************************** C=* From: SPIDER - MODULAR IMAGE PROCESSING SYSTEM * C=* Copyright (C)2002,2009 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 C PURPOSE: FOR: 'BPD 32F' IMPROVED WITH LESS MEMORY REQUIRED AND C MORE STACK USAGE. 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/ CALL SET_MPI(ICOMM,MYPID,MPIERR) ! SETS ICOMM AND MYPID 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 REVERSE 3D FFT ON XE USING ALL THREADS c write(6,*) ' wiw32d; inv plan on xe: ',n2 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 c write(6,*) 'wiw32d; nrmw2 xe : ',nsam,n2 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 C write(6,*) ' wiw32d; nrmw2 xe : ',nsam,n2 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 C write(6,*) ' wiw32d; nrmw2 xo : ',nsam,n2 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 CHARACTER(LEN=*) :: FILPAT 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,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. INCLUDE 'mpif.h' INTEGER :: ISTAT(MPI_STATUS_SIZE) INTEGER, ALLOCATABLE, DIMENSION(:) :: PSIZE, NBASE REAL , ALLOCATABLE, DIMENSION(:,:,:) :: PRJLOC REAL , ALLOCATABLE, DIMENSION(:,:,:) :: PRJBUF REAL , ALLOCATABLE, DIMENSION(:,:,:) :: WELOC, WOLOC COMPLEX, ALLOCATABLE, DIMENSION(:,:,:) :: XELOC, XOLOC ICOMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(ICOMM, MYPID, MPIERR) CALL MPI_COMM_SIZE(ICOMM, 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 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. C CREATE FORWARD FFTW3 PLAN FOR 2D FFT ON BI USING ONE THREAD c write(6,*) ' wiw32d; inv plan on bi,n : ',n CALL FMRS_PLAN(.TRUE.,BI,N,N,1, 1,+1,IRTFLG) IF (IRTFLG .NE. 0) RETURN 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. C PARTITIONS NANG DATA ITEMS EVENLY AMONG AVAILABLE PROCESSORS. C PSIZE HAS NUMBER OF PARTICLES, NBASE HAS STARTING LOCATION. CALL SETPART(NANG, PSIZE, NBASE) C NANGLOC IS THE # OF PARTICLES TO BE USED FOR CURRENT PROCESS. 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 SO EACH PROCESSOR GETS DIFFERENT SET OF IMAGE FILES C (ASSUMES SUFFICIENT MEMORY TO HOLD NANG/NPROCS IMAGES) C FOR IPROC =1, MYPID=0 READS AND KEEPS IMAGES FOR ITSELF C FOR IPROC =1, MYPID>0 SKIPS READ C FOR IPROC >1, MYPID=0 READS AND BROADCASTS IMAGES C FOR IPROC >1, MYPID>0 SKIPS READ AND RECEIVES IMAGES DO IPROC = 1, NPROCS NANGLOC = PSIZE(IPROC) C READ IMAGES INTO THE PRJBUF BUFFER FIRST, THEN DISTRIBUTE DO JLOC = 1, NANGLOC JGLB = NBASE(IPROC) + JLOC NLET = 0 ! TO USE: lnblnkn IN FILGET CALL FILGET(FILPAT,FILNAM,NLET,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 C ONLY PROCESS WITH MYPID = 0 NEEDS TO READ IF (MYPID .LE. 0) THEN DO IROW = 1,NROW CALL REDLIN1P(INPROJ,PRJBUF(1,IROW,JLOC),NSAM,IROW) ENDDO ENDIF CLOSE(INPROJ) ENDDO IF (IPROC .GT. 1) THEN #ifdef MPI_DEBUG write(6,*) 'wiw32dq: iproc = ', iproc #endif IF (MYPID .EQ. 0) THEN C SEND PROJECTION IMAGES TO OTHER PROCESSOR(S) #ifdef MPI_DEBUG write(6,*) 'wiw32dq: sending to mypid= ', mypid #endif CALL MPI_SEND(PRJBUF , NSAM*NSAM*NANGLOC, MPI_REAL, & IPROC-1, IPROC-1, ICOMM, 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 #ifdef MPI_DEBUG WRITE(6,*) 'wiw32dq: receiving mypid: ', mypid #endif CALL MPI_RECV(PRJLOC, NSAM*NSAM*NANGLOC, MPI_REAL, & 0, MPI_ANY_TAG, ICOMM, 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 PROJECTION IMAGES FOR THIS PROCESS #ifdef MPI_DEBUG write(6,*) 'wiw32dq: keeping mypid= ', mypid #endif 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 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 c write(6,*) ' wiw32dq: xeloc zero = ', mypid 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 INV = +1 ! FORWARD FOURIER TRANSFORMS 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 FORWARD FOURIER TRANSFORM OF: BI c write(6,*) ' wiw32d; fft on bi,n : ',n 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 c write(6,*) ' before reduce ; n,nsam,nangloc: ', n,nsam,nangloc c write(6,*) ' xe:' , xeloc(0,1,1),xeloc(0,1,2),xeloc(75,75,75) c write(6,*) ' xo:' , xoloc(0,1,1),xoloc(0,1,2),xoloc(75,75,75) c write(6,*) ' we:' , weloc(0,1,1),weloc(0,1,2),weloc(75,75,75) c write(6,*) ' wo:' , woloc(0,1,1),woloc(0,1,2),woloc(75,75,75) C SUM UP X AND W FROM THE LOCAL PIECES (XLOC, WLOC) C RESIDING ON EACH PROCESSOR NVAL = (NSAM + 1) * N DO K3 = 1, N CALL ALLREDUCE_MPI('WIW32DQ','XELOC', & XELOC(0,1,K3),XE(0,1,K3),NVAL, 'X','S',ICOMM) CALL ALLREDUCE_MPI('WIW32DQ','XOLOC', & XOLOC(0,1,K3),XO(0,1,K3),NVAL, 'X','S',ICOMM) CALL ALLREDUCE_MPI('WIW32DQ','WELOC', & WELOC(0,1,K3),WE(0,1,K3),NVAL, 'R','S',ICOMM) CALL ALLREDUCE_MPI('WIW32DQ','WOLOC', & WOLOC(0,1,K3),WO(0,1,K3),NVAL, 'R','S',ICOMM) ENDDO IF (ALLOCATED(PRJLOC)) DEALLOCATE(PRJLOC) 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 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 #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 FORWARD 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, 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 #endif C -------------------- END OF NON-MPI VERSION ------------------ c write(6,*) ' ' c write(6,*) ' xe:' , xe(0,1,1),xe(0,1,2),xe(75,75,75) c write(6,*) ' xo:' , xo(0,1,1),xo(0,1,2),xo(75,75,75) c write(6,*) ' we:' , we(0,1,1),we(0,1,2),we(75,75,75) c write(6,*) ' wo:' , wo(0,1,1),wo(0,1,2),wo(75,75,75) 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 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