C++********************************************************************* C C REPCG.F CORRECTIONS APPLIED ON VOLUME SIDE 11/09/98 C COMPRESSION OF ANGLES 08/14/96 C USED PROJT FOR BCKCQ CALL FEB 2000 ARDEAN LEITH C OPFILEC FEB 03 ARDEAN LEITH C VERBOSE FEB 06 ARDEAN LEITH C MPI BUG FIXED OCT 08 ARDEAN LEITH C REFACTORED OCT 08 ARDEAN LEITH C ********************************************************************** C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM. AUTHOR: J.FRANK * C=* Copyright (C) 1985-2008 Health Research Inc. * C=* * C=* HEALTH RESEARCH INCORPORATED (HRI), * C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455. * C=* * C=* Email: spider@wadsworth.org * C=* * C=* This program is free software; you can redistribute it and/or * C=* modify it under the terms of the GNU General Public License as * C=* published by the Free Software Foundation; either version 2 of the * C=* License, or (at your option) any later version. * C=* * C=* This program is distributed in the hope that it will be useful, * C=* but WITHOUT ANY WARRANTY; without even the implied warranty of * C=* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * C=* General Public License for more details. * C=* * C=* You should have received a copy of the GNU General Public License * C=* along with this program; if not, write to the * C=* Free Software Foundation, Inc., * C=* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * C=* * C ********************************************************************** C C REPCG C C PURPOSE: CALCULATES 3D RECONSTRUCTION USING CONJUGATE GRADIENTS C WITH REGULARIZATION. C THE RECONSTRUCTION IS CALCULATED INSIDE A SPHERE OF C SPECIFIED RADIUS C RECONSTRUCTION KEPT IN SQUARE TO INTRODUCE OTHER C CONSTRAINTS. AVERAGE OUTSIDE WINDOW IS SUBTRACTED C SYMMETRIES NOT IMPLEMENTED C C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* SUBROUTINE REPCG INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' INCLUDE 'F90ALLOC.INC' REAL, DIMENSION(:,:), POINTER :: PANG REAL, ALLOCATABLE, DIMENSION(:,:,:) :: CB REAL, ALLOCATABLE, DIMENSION(:,:) :: ANG,DM REAL, ALLOCATABLE, DIMENSION(:) :: BCKN INTEGER, ALLOCATABLE, DIMENSION(:,:) :: QM INTEGER, ALLOCATABLE, DIMENSION(:) :: ILIST,LB LOGICAL :: MD CHARACTER(LEN=MAXNAM) :: ANGDOC,FILNAM,FINPAT CHARACTER(LEN=1) :: NULL,ANS COMMON /PAR/ LDP,NM,LDPNM DATA INPIC/99/ DATA INPROJ/98/,LUNUNK/77/ #ifdef USE_MPI include 'mpif.h' ICOMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(ICOMM, MYPID, MPIERR) #else MYPID = -1 #endif NULL = CHAR(0) NILMAX = NIMAX C N - LINEAR DIMENSION OF PROJECTIONS AND RESTORED CUBE C NANG - NUMBER OF ANGLES (PROJECTIONS) ALLOCATE (ILIST(NILMAX), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'REPCG, ILIST',NILMAX) RETURN ENDIF CALL FILELIST(.TRUE.,INPIC,FINPAT,NLET,ILIST,NILMAX,NANG, & 'ENTER TEMPLATE FOR 2-D IMAGES',IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(27,'REPCG ',NE) GOTO 9999 ENDIF MAXNUM = MAXVAL(ILIST(1:NANG)) C NANG - TOTAL NUMBER OF IMAGES IF (MYPID .LE. 0) WRITE(NOUT,2001) NANG 2001 FORMAT(' NUMBER OF IMAGES=',I6) CALL RDPRM(RI,NOT_USED,'RADIUS OF RECONSTRUCTED OBJECT') C RETRIEVE ARRAY WITH ANGLES DATA IN IT MAXXT = 4 MAXYT = MAXNUM CALL GETDOCDAT('ANGLES DOC',.TRUE.,ANGDOC,LUNUNK,.FALSE.,MAXXT, & MAXYT,PANG,IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(4,'REPCG ',NE) GOTO 9999 ENDIF C OPEN FIRST IMAGE FILE TO DETERMINE NSAM, NROW, NSL CALL FILGET(FINPAT,FILNAM,NLET,ILIST(1),IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(4,'REPCG ',NE) GOTO 9999 ENDIF MAXIM = 0 CALL OPFILEC(0,.FALSE.,FILNAM,INPIC,'O',IFORM,NSAM,NROW,NSL, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 N = NSAM M = N NM = (N-M) / 2 LDP = M / 2 + 1 LDPNM = LDP + NM CALL RDPRMC (ANS, NLETI, .TRUE., & 'DOES YOUR VOLUME HAVE SYMMETRIES? (NOT IMPLEMENTED) (Y/N)', & NULL,IRT) IF (ANS .EQ. 'Y') THEN CALL ERRT(101,'*** SYMMETRIES OPTION NOT IMPLEMENTED',NE) GOTO 9999 ENDIF C DUM IS A DUMMY VARIABLE, VALUE OF NN IS DETERMINED MD = .FALSE. CALL PREPCUB_S(N,NN,DUM,RI,MD) C USE NN TO ALLOCATE (QM) ALLOCATE (QM(5,NN), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'REPS, QM',5*NN) GOTO 9999 ENDIF MD = .TRUE. CALL PREPCUB_S(N,NN,QM,RI,MD) NMAT = N*N*N LSD = N+2-MOD(N,2) C IN THIS VERSION TOTAL MEMORY IS THREE VOLUMES BCKN, BCKE AND CUBE C PLUS 2 2D PROJECTIONS. C CUBE - KEEPS BACK-PROJECTED ORIGINAL PROJECTIONS, READ FROM THE DISK C BCKE - WORKING VOLUME C BCKN - CURRENT RECONSTRUCTION C FIND NUMBER OF OMP THREADS CALL GETTHREADS(NUMTH) MAXIM = 0 IFORM = 3 CALL OPFILEC(0,.TRUE.,FILNAM,INPIC,'U',IFORM,M,M,M, & MAXIM,'RECONSTRUCTED 3-D',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ALLOCATE(CB(N,N,N), & ANG(3,NANG), & DM(9,NANG), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'REPCG, CB, ANG & DM',IER) GOTO 9999 ENDIF CALL REDPRCG(N,NANG,CB,ANG,ILIST,QM,NN,DM,RI,PANG, & FINPAT,BNORM,CHI2) C COMPRESS ANGLES - IT CHANGES NANG !! ALLOCATE (LB(NANG), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'REPCG, LB',IER) GOTO 9999 ENDIF C CONVERTS ANGLES TO 'DM' FORMAT? CALL HIANG(ANG,NANG,DM,LB,LO) NANG = LO IF (MYPID .LE. 0) WRITE(NOUT,2027) NANG 2027 FORMAT(' EFFECTIVE NUMBER OF ANGLES: ',I6) ALLOCATE(BCKN(N*N*N), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'REPCG, BCKN',N*N*N) GOTO 9999 ENDIF CALL REPRCG3_Q(CB,BCKN,LSD,N,QM,NN,DM,LB,NANG,RI, & NUMTH,BNORM,CHI2) DO K=1,M*M MI = (K-1)*M CALL WRTLIN(INPIC,BCKN(1+MI),M,K) ENDDO 9999 IF (ALLOCATED(BCKN)) DEALLOCATE(BCKN) IF (ALLOCATED(LB)) DEALLOCATE(LB) IF (ALLOCATED(DM)) DEALLOCATE(DM) IF (ALLOCATED(ANG)) DEALLOCATE(ANG) IF (ALLOCATED(CB)) DEALLOCATE(CB) IF (ALLOCATED(QM)) DEALLOCATE(QM) IF (ASSOCIATED(PANG)) DEALLOCATE(PANG) IF (ALLOCATED(ILIST))DEALLOCATE(ILIST) CLOSE(INPIC) END C++********************************************************************* C--********************************************************************* SUBROUTINE REDPRCG(N,NANG,CB,ANG,ILIST,IPCUBE,NN,DM,RI,ANGBUF, & FINPAT,BNORM,CHI2) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' REAL :: CB(N,N,N),PROJ(N,N),ANG(3,NANG) INTEGER :: ILIST(NANG),IPCUBE(5,NN) REAL :: DM(9,NANG) REAL :: ANGBUF(4,NANG) DOUBLE PRECISION :: ABA,SUS,SSQ CHARACTER(LEN=MAXNAM) :: FILNAM,FINPAT DATA INPROJ/98/ #ifdef USE_MPI include 'mpif.h' REAL, ALLOCATABLE, DIMENSION(:,:,:) :: CB_LOC, PRJBUF, PRJLOC REAL, ALLOCATABLE, DIMENSION(:,:) :: ANG_LOC INTEGER, ALLOCATABLE, DIMENSION(:) :: PSIZE INTEGER, ALLOCATABLE, DIMENSION(:) :: NBASE DOUBLE PRECISION :: ABA_LOC,SUS_LOC,SSQ_LOC INTEGER :: MPISTAT(MPI_STATUS_SIZE) ICOMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(ICOMM, MYPID , MPIERR) CALL MPI_COMM_SIZE(ICOMM, NPROCS, MPIERR) ALLOCATE(PSIZE(NPROCS),NBASE(NPROCS), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,' REPRCG3_Q, PSIZE & NBASE.',2*NPROCS) STOP ENDIF #else MYPID = -1 #endif ABA = 0.0D0 KLP = 0 SUS = 0.0D0 SSQ = 0.0D0 KLS = 0 c$omp parallel do private(i,j,k) DO K=1,N DO J=1,N DO I=1,N CB(I,J,K) = 0.0 ENDDO ENDDO ENDDO #ifdef USE_MPI ABA_LOC = 0.0D0 KLP_LOC = 0 KLS_LOC = 0 SUS_LOC = 0.0D0 SSQ_LOC = 0.0D0 CALL SETPART(NANG, PSIZE, NBASE) NANGLOC = PSIZE(MYPID+1) ALLOCATE(PRJBUF(N,N,PSIZE(1)), & PRJLOC(N,N,NANGLOC), & ANG_LOC(3,NANGLOC), & CB_LOC(N,N,N), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = N*N*PSIZE(1) + N*N*NANGLOC + 3*NANGLOC + N*N*N CALL ERRT(46,'REPRCG3_Q, PRJBUF, PRJLOC..',MWANT) RETURN ENDIF C ZERO LOCAL ARRAYS CB_LOC = 0.0 ANG_LOC = 0.0 DO IPROC = 1, NPROCS NLOC = PSIZE(IPROC) C READ A SUBSET OF IMAGES (ONLY ONE PROCESSOR READS) DO K=1,NLOC KGLB = K + NBASE(IPROC) NLET = 0 CALL FILGET(FINPAT,FILNAM,NLET,ILIST(KGLB),IRTFLG) IF (IRTFLG .NE. 0) RETURN MAXIM = 0 CALL OPFILEC(0,.FALSE.,FILNAM,INPROJ,'O',IFORM, & LSAM,LROW,NSL, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN C DO K2=1,N CALL REDLIN1P(INPROJ,PRJBUF(1,K2,K),N,K2) ENDDO IF (MYPID .EQ. 0) CLOSE(INPROJ) ENDDO C DISTRIBUTE IMAGES IF (IPROC .GT. 1) THEN IF (MYPID .EQ. 0) THEN CALL SEND_MPI('REDPRCG','PRJBUF', PRJBUF, N*N*NLOC, & 'R',IPROC-1,IPROC-1, ICOMM) ELSE IF (MYPID .EQ. IPROC-1) THEN CALL MPI_RECV(PRJLOC, N*N*NLOC, MPI_REAL, & 0, MPI_ANY_TAG, ICOMM, & MPISTAT, MPIERR) IF (MPIERR .NE. 0) THEN WRITE(6,*) 'REDPRCG: RECV FAILED' STOP ENDIF ENDIF ELSE IF (MYPID .EQ. 0) THEN CALL SCOPY(N*N*NLOC,PRJBUF,1,PRJLOC,1) ENDIF ENDDO DO K = 1, NANGLOC KGLB = K + NBASE(MYPID+1) C ORDER IN DOCUMENT FILE IS PSI, THETA, PHI AND ANGLES ARE IN C DEGREES! C IN ANG ARRAY IT IS OTHER WAY AROUND ITMP = ILIST(KGLB) ANG_LOC(3,K) = ANGBUF(2,ITMP) ANG_LOC(2,K) = ANGBUF(3,ITMP) ANG_LOC(1,K) = ANGBUF(4,ITMP) c write(6,91) kglb,itmp,(ang_loc(j,k),j=3,1,-1) 91 format(' kglb:',i5,' proj:',i6, & '; psi:',f6.1,' theta:',f6.1,' phi:',f6.1) C ESTIMATE AVERAGE OUTSIDE THE CIRCLE? CALL ASTASQ(PRJLOC(1,1,K), N, RI, ABA_LOC, KLP_LOC, & SUS_LOC, SSQ_LOC, KLS_LOC) C BACKPROJECTION? CALL RPRQ(N,PRJLOC(1,1,K),CB_LOC,IPCUBE,NN, & ANG_LOC(1,K),ANG_LOC(2,K),ANG_LOC(3,K),DM(1,KGLB),RI) ENDDO IF (ALLOCATED(PRJBUF)) DEALLOCATE(PRJBUF) IF (ALLOCATED(PRJLOC)) DEALLOCATE(PRJLOC) C GATHER ANG_LOC INTO ANG DO J = 1, NPROCS NBASE(J) = 3 * NBASE(J) PSIZE(J) = 3 * PSIZE(J) ENDDO CALL MPI_ALLGATHERV(ANG_LOC, 3*NANGLOC, MPI_REAL, & ANG, PSIZE, NBASE, & MPI_REAL, ICOMM, IERR) DO K = 1, NANG IF (VERBOSE) THEN IF (MYPID .EQ. 0) WRITE(NOUT,3331) K,(ANG(J,K),J=3,1,-1) 3331 FORMAT(' PROJECTION #',I6, & '; PSI=',F6.1,' THETA=',F6.1,' PHI=',F6.1) ENDIF END DO N3 = N*N*N IERR = 0 CALL ALLREDUCE_MPI('REPRCG3_Q','CB', CB_LOC,CB, & N3, 'R','S',ICOMM) CALL ALLREDUCE_MPI('REPRCG3_Q','ABA', ABA_LOC,ABA, & 1, 'D','S',ICOMM) CALL ALLREDUCE_MPI('REPRCG3_Q','KLP', KLP_LOC,KLP, & 1, 'I','S',ICOMM) CALL ALLREDUCE_MPI('REPRCG3_Q','SSQ', SSQ_LOC,SSQ, & 1, 'D','S',ICOMM) CALL ALLREDUCE_MPI('REPRCG3_Q','SUS', SUS_LOC,SUS, & 1, 'D','S',ICOMM) CALL ALLREDUCE_MPI('REPRCG3_Q','KLS', KLS_LOC,KLS, & 1, 'I','S',ICOMM) IF (ALLOCATED(CB_LOC)) DEALLOCATE(CB_LOC) IF (ALLOCATED(ANG_LOC)) DEALLOCATE(ANG_LOC) IF (ALLOCATED(PSIZE)) DEALLOCATE(PSIZE) IF (ALLOCATED(NBASE)) DEALLOCATE(NBASE) C --------------- END OF MPI CODE ------------------------------ #else DO K=1,NANG NLET = 0 CALL FILGET(FINPAT,FILNAM,NLET,ILIST(K),IRTFLG) IF (IRTFLG .NE. 0) RETURN MAXIM = 0 CALL OPFILEC(0,.FALSE.,FILNAM,INPROJ,'O',IFORM, & LSAM,LROW,NSL, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN DO K2=1,N CALL REDLIN(INPROJ,PROJ(1,K2),N,K2) ENDDO CLOSE(INPROJ) C ORDER IN DOCUMENT FILE IS PSI, THETA, PHI AND ANGLES ARE IN C DEGREES! IN ANG ARRAY IT IS OTHER WAY AROUND ITMP = ILIST(K) ICOUNT = ANGBUF(1,ITMP) IF (ICOUNT .LE. 0) THEN C MISSING KEY CALL ERRT(102,'MISSING ANGLE FOR IMAGE',ITMP) RETURN ENDIF ANG(3,K) = ANGBUF(2,ITMP) ANG(2,K) = ANGBUF(3,ITMP) ANG(1,K) = ANGBUF(4,ITMP) IF (VERBOSE) THEN WRITE(NOUT,3331) K,(ANG(J,K),J=3,1,-1) 3331 FORMAT(' PROJECTION #',I6, & '; PSI=',F6.1,' THETA=',F6.1,' PHI=',F6.1) ENDIF C ESTIMATE AVERAGE OUTSIDE THE CIRCLE? CALL ASTASQ(PROJ,N,RI,ABA,KLP,SUS,SSQ,KLS) C BACKPROJECTION? CALL RPRQ(N,PROJ,CB,IPCUBE,NN, & ANG(1,K),ANG(2,K),ANG(3,K),DM(1,K),RI) ENDDO #endif C --------------- END OF NON-MPI CODE ------------------------- ABA = ABA / KLP C PRINT STATISTICS IF (MYPID .LE. 0) WRITE(NOUT,2044) KLP,ABA 2044 FORMAT & (' TOTAL NUMBER OF POINTS IN PROJECTIONS =',I10, & /,' AVERAGE OUTSIDE THE WINDOW =',1PE10.3,/) C SUBTRACT AVERAGE AND WRITE TO THE FILE BNORM = 0.0 QT = ABA * NANG c$omp parallel do private(kn,j,k,i),reduction(+:bnorm) DO KN=1,NN J = IPCUBE(4,KN) K = IPCUBE(5,KN) DO I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN) CB(I,J,K) = CB(I,J,K) - QT BNORM = BNORM + CB(I,J,K) * CB(I,J,K) ENDDO ENDDO C CALCULATE SUM OF SQUARES OF PROJECTIONS WITHIN CIRCLES AFTER C SUBTRACTION OF THE AVERAGE SSQ = SSQ - ABA * (2 * SUS - ABA * KLS) CHI2 = SSQ IF (MYPID .LE. 0) WRITE(NOUT,2045) BNORM,CHI2 2045 FORMAT (' SQUARED BP =',1PE10.3, & /,' CHI2 =',1PE10.3,/) END C++********************************************************************* C--********************************************************************* SUBROUTINE REPRCG3_Q(CB,BCKN,LSD,N,IPCUBE,NN,DM,LB,NANG,RI, & NUMTH,BNORM,CHI2) C NUMTH - NUMTHREDS FOR MP; =NUMTHREADS() FOR MP, OTHERWISE=1. INCLUDE 'CMBLOCK.INC' REAL :: CB(N,N,N),BCKN(N,N,N) REAL :: DM(9,NANG),PROJ(N*N,NUMTH) REAL :: PROJT(4,N*N) INTEGER :: IPCUBE(5,NN),LB(NANG) REAL,ALLOCATABLE, DIMENSION(:,:,:) :: BCKE,BCKP COMMON /PAR/ LDP,NM,LDPNM #ifdef USE_MPI include 'mpif.h' INTEGER :: MPISTAT(MPI_STATUS_SIZE) REAL :: ALA, SQ, DELSQ INTEGER,ALLOCATABLE, DIMENSION(:) :: PSIZE, NBASE INTEGER,ALLOCATABLE, DIMENSION(:) :: LBLOC REAL, ALLOCATABLE, DIMENSION(:,:) :: DMLOC REAL, ALLOCATABLE, DIMENSION(:,:,:):: BCKE_SUM REAL, ALLOCATABLE, DIMENSION(:,:,:):: RESID ICOMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(ICOMM, MYPID , MPIERR) CALL MPI_COMM_SIZE(ICOMM, NPROCS, MPIERR) ALLOCATE(PSIZE(NPROCS), NBASE(NPROCS), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = 2 * NPROCS CALL ERRT(46,'REPRCG3_Q, NBASE',MWANT) RETURN ENDIF CALL SETPART(NANG, PSIZE, NBASE) NANGLOC = PSIZE(MYPID+1) #else MYPID = -1 #endif CALL RDPRM2(ERRM,CHIM,NOT_USED,'ERROR LIMIT, CHI^2 LIMIT') CALL RDPRMI(MAXIT,MODE,NOT_USED,'ITERATION LIMIT, MODE') CALL RDPRM(ALA,NOT_USED,'LAMBDA') DI = N / 2 IF (MODE .EQ. 1 .AND. RI+1.0 .GT. DI & .OR. MODE .EQ. 2 .AND. RI+2.0 .GT. DI & .OR. MODE .EQ. 3 .AND. RI+3.0 .GT. DI) THEN CALL ERRT(101,'RADIUS TOO LARGE',NE) GOTO 9999 ENDIF ALLOCATE (BCKE(N,N,N),BCKP(N,N,N), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN MWANT = 2*N*N*N CALL ERRT(46,'REPRCG3_Q, BCKE,BCKP',MWANT) GOTO 9999 ENDIF #ifdef USE_MPI ALLOCATE(BCKE_SUM(N,N,N), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN MWANT = N*N*N CALL ERRT(46,'REPRCG3_Q, BCKE_SUM',MWANT) GOTO 9999 ENDIF #endif NMAT = N*N*N LTB = N*N c$omp parallel do private(k,j,i) DO K=1,N DO J=1,N DO I=1,N BCKN(I,J,K) = 0.0 BCKP(I,J,K) = 0.0 #ifdef USE_MPI BCKE_SUM(I,J,K) = 0.0 #endif ENDDO ENDDO ENDDO ERR = 1.0 ITER = 0 Q = 0.0 C ------------------ START OF: MPI CODE ------------------------ #ifdef USE_MPI C DISTRIBUTE LB AND DM ALLOCATE(LBLOC(NANGLOC), & DMLOC(9, NANGLOC), & STAT = IRTFLAG) IF (IRTFLG .NE. 0) THEN MWANT = 10 * NANGLOC CALL ERRT(46,'REPRCG3_Q, DMLOC',MWANT) GOTO 9999 ENDIF IDMTAG = 1 LBTAG = 2 MASTER = 0 IF (MYPID .EQ. 0) THEN C PROCESS 0 DISTRIBUTES DM() AND LB() TO OTHER PROCESSORS DO IP = 2, NPROCS IBEGIN = NBASE(IP) + 1 CALL SEND_MPI('REPRCG3_Q','DM', DM(1,IBEGIN), 9*NANGLOC, & 'R',IP-1,IDMTAG,ICOMM) CALL SEND_MPI('REPRCG3_Q','LB', LB(IBEGIN), NANGLOC, & 'I',IP-1,LBTAG,ICOMM) ENDDO CALL SCOPY(9*NANGLOC, DM, 1, DMLOC, 1) CALL SCOPY(NANGLOC , LB, 1, LBLOC, 1) ELSE C === SLAVES RECEIVE DATA FROM THE MASTER === CALL MPI_RECV(DMLOC, 9*NANGLOC, MPI_REAL, MASTER, & IDMTAG, ICOMM, MPISTAT, IERR) CALL MPI_RECV(LBLOC, NANGLOC, MPI_INTEGER, MASTER, & LBTAG, ICOMM, MPISTAT, IERR) ENDIF #ifdef MPI_DEBUG c write(6,*)' reprcg3_q: data distributed, mypid:',mypid #endif DO ITER=1,MAXIT DELSQ = 0.0 c$omp parallel do private(kn,j,k,i),reduction(+:delsq) DO KN=1,NN J = IPCUBE(4,KN) K = IPCUBE(5,KN) DO I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN) DELSQ = DELSQ + CB(I,J,K) * CB(I,J,K) ENDDO ENDDO Q = Q * DELSQ c$omp parallel do private(kn,j,k,i) DO KN=1,NN J = IPCUBE(4,KN) K = IPCUBE(5,KN) DO I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN) BCKP(I,J,K) = CB(I,J,K) - Q * BCKP(I,J,K) ENDDO ENDDO Q = -1.0 / DELSQ c$omp parallel do private(k,j,i) DO K=1,N DO J=1,N DO I=1,N BCKE(I,J,K) = 0.0 BCKE_SUM(I,J,K) = 0.0 ENDDO ENDDO ENDDO C BCKP -> PROJ -> BCKE C LOOP OVER PROJECTIONS DO K=1,NANGLOC L_EN = MIN(NANGLOC,K+NUMTH-1) L_NUM = MIN(NUMTH,NANGLOC-K+1) c$omp parallel do private(i,j) DO J=1,L_NUM DO I=1,LTB PROJ(I,J) = 0.0 ENDDO ENDDO C PROJECT CALL PRJCQ(BCKP,NMAT,DMLOC(1,K), & PROJ,N,IPCUBE,NN) C HERE BCKCQ ITSELF IS MP DO L_TH=1,L_NUM C MULTIPLY PROJECTIONS BY THEIR WEIGHTS IF (LBLOC(K+L_TH-1) .GT. 1) THEN c$omp parallel do private(i) DO I=1,N*N PROJ(I,L_TH) = PROJ(I,L_TH) * LBLOC(K+L_TH-1) ENDDO ENDIF c$omp parallel do private(i,pt) DO I = 1,N*N - N - 1 PT = PROJ(I, L_TH) PROJT(1,I) = PT PROJT(2,I) = PROJ(I+N, L_TH) - PT PROJT(3,I) = PROJ(I+1, L_TH) - PT PROJT(4,I) = PROJ(I+N+1,L_TH) - PROJ(I+1,L_TH) - & PROJT(2,I) ENDDO C BACKPROJECT CALL BCKCQ(BCKE,NMAT,DMLOC(1,K+L_TH-1), & PROJT,N,IPCUBE,NN) ENDDO ENDDO C CALL ALLREDUCE_MPI('REPRCG3_Q','BCKE_SUM',BCKE,BCKE_SUM, & NMAT, 'R','S',ICOMM) C === IGNORE MODE > 0 FOR THE TIME BEING === IF (MODE .EQ. 1) THEN CALL FIXEDGE1(BCKP,NMAT,BCKP,N,IPCUBE,NN) c CALL BFIRSTS(BCKE,BCKP,N,ALA,IPCUBE,NN,RI) sep 01 al CALL BFIRSTS(BCKE_SUM,BCKP,N,ALA,IPCUBE,NN) ELSEIF(MODE.EQ.2) THEN CALL FIXEDGE2(BCKP,NMAT,BCKP,N,IPCUBE,NN) CALL BSECOND(BCKE_SUM,BCKP,N,ALA,IPCUBE,NN) ELSEIF(MODE.EQ.3) THEN CALL FIXEDGE3(BCKP,NMAT,BCKP,N,IPCUBE,NN) CALL BTHIRD(BCKE_SUM,BCKP,N,ALA,IPCUBE,NN) ENDIF AKDEN = 0.0 c$omp parallel do private(kn,j,k,i),reduction(+:akden) DO KN=1,NN J = IPCUBE(4,KN) K = IPCUBE(5,KN) DO I=IPCUBE(3,KN), & IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN) AKDEN = AKDEN+BCKP(I,J,K) * BCKE_SUM(I,J,K) ENDDO ENDDO P=DELSQ / AKDEN c$omp parallel do private(kn,j,k,i) DO KN=1,NN J = IPCUBE(4,KN) K = IPCUBE(5,KN) DO I=IPCUBE(3,KN), & IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN) BCKN(I,J,K) = BCKN(I,J,K)+P*BCKP(I,J,K) CB(I,J,K) = CB(I,J,K)-P*BCKE_sum(I,J,K) ENDDO ENDDO ERR = DELSQ / BNORM CHI2 = CHI2 - P * DELSQ IF (MYPID .EQ. 0) WRITE(NOUT,2041) ITER,ERR,CHI2 2041 FORMAT(' ITERATION:', I3,' ERROR:',1PE12.4, & ' CHISQ:',1PE12.4) C CHECK STOPPING CRITERIA IF (ABS(ERR) .LE. ERRM .OR. ABS(CHI2) .LE. CHIM) GOTO 9999 ENDDO C END OF: -------------- USE_MPI ------------------------------ #else DO ITER=1,MAXIT DELSQ = 0.0 c$omp parallel do private(kn,j,k,i),reduction(+:DELSQ) DO KN=1,NN J = IPCUBE(4,KN) K = IPCUBE(5,KN) DO I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN) DELSQ=DELSQ+CB(I,J,K)*CB(I,J,K) ENDDO ENDDO Q = Q * DELSQ c$omp parallel do private(kn,j,k,i) DO KN=1,NN J = IPCUBE(4,KN) K = IPCUBE(5,KN) DO I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN) BCKP(I,J,K) = CB(I,J,K) - Q * BCKP(I,J,K) ENDDO ENDDO Q = -1.0 / DELSQ c$omp parallel do private(k,j,i) DO K=1,N DO J=1,N DO I=1,N BCKE(I,J,K) = 0.0 ENDDO ENDDO ENDDO C BCKP -> PROJ -> BCKE C LOOP OVER PROJECTIONS DO K=1,NANG,NUMTH L_EN = MIN(NANG,K+NUMTH-1) L_NUM = MIN(NUMTH,NANG-K+1) c$omp parallel do private(i,j) DO J=1,L_NUM DO I=1,LTB PROJ(I,J) = 0.0 ENDDO ENDDO c$omp parallel do private(l_th),shared(nmat,n,nn),schedule(static) DO L_TH=1,L_NUM CALL PRJCQ(BCKP,NMAT,DM(1,K+L_TH-1), & PROJ(1,L_TH),N,IPCUBE,NN) ENDDO C HERE BCKCQ ITSELF IS MP DO L_TH=1,L_NUM C MULTIPLY PROJECTIONS BY THEIR WEIGHTS IF (LB(K+L_TH-1) .GT. 1) THEN c$omp parallel do private(i) DO I=1,N*N PROJ(I,L_TH) = PROJ(I,L_TH)*LB(K+L_TH-1) ENDDO ENDIF c$omp parallel do private(i,pt) DO I = 1,N*N - N - 1 PT = PROJ(I, L_TH) PROJT(1,I) = PT PROJT(2,I) = PROJ(I+N, L_TH) - PT PROJT(3,I) = PROJ(I+1, L_TH) - PT PROJT(4,I) = PROJ(I+N+1,L_TH) - PROJ(I+1,L_TH) - & PROJT(2,I) ENDDO C BACKPROJECT CALL BCKCQ(BCKE,NMAT,DM(1,K+L_TH-1),PROJT,N,IPCUBE,NN) ENDDO ENDDO IF (MODE .EQ. 1) THEN CALL FIXEDGE1(BCKP,NMAT,BCKP,N,IPCUBE,NN) c CALL BFIRSTS(BCKE,BCKP,N,ALA,IPCUBE,NN,RI) sep 01 al CALL BFIRSTS(BCKE,BCKP,N,ALA,IPCUBE,NN) ELSEIF(MODE .EQ. 2) THEN CALL FIXEDGE2(BCKP,NMAT,BCKP,N,IPCUBE,NN) CALL BSECOND(BCKE,BCKP,N,ALA,IPCUBE,NN) ELSEIF (MODE .EQ. 3) THEN CALL FIXEDGE3(BCKP,NMAT,BCKP,N,IPCUBE,NN) CALL BTHIRD(BCKE,BCKP,N,ALA,IPCUBE,NN) ENDIF AKDEN = 0.0 c$omp parallel do private(kn,j,k,i),reduction(+:akden) DO KN=1,NN J = IPCUBE(4,KN) K = IPCUBE(5,KN) DO I=IPCUBE(3,KN), & IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN) AKDEN = AKDEN + BCKP(I,J,K) * BCKE(I,J,K) ENDDO ENDDO P = DELSQ / AKDEN c$omp parallel do private(kn,j,k,i) DO KN=1,NN J = IPCUBE(4,KN) K = IPCUBE(5,KN) DO I=IPCUBE(3,KN), & IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN) BCKN(I,J,K) = BCKN(I,J,K)+P*BCKP(I,J,K) CB(I,J,K) = CB(I,J,K)-P*BCKE(I,J,K) ENDDO ENDDO ERR = DELSQ / BNORM CHI2 = CHI2 - P * DELSQ IF (VERBOSE) THEN C SKIP THIS IF NOT VERBOSE WRITE(NOUT,2041) ITER,ERR,CHI2 2041 FORMAT(' ITERATION: 'I3,' ERROR:',1PE12.4,' & CHISQ: ',1PE12.4) ENDIF C CHECK STOPPING CRITERIA IF (ABS(ERR) .LE. ERRM .OR. ABS(CHI2) .LE. CHIM) THEN C WRITE FINAL ERROR IF (.NOT. VERBOSE) WRITE(NOUT,2041) ITER,ERR,CHI2 GOTO 9999 ENDIF ENDDO C END OF: -------------- NON-MPI ---------------------------- #endif C NUMBER OF ITERATIONS EXCEEDED 9999 IF (ALLOCATED(BCKE)) DEALLOCATE (BCKE) IF (ALLOCATED(BCKP)) DEALLOCATE (BCKP) #ifdef USE_MPI IF (ALLOCATED(BCKE_SUM)) DEALLOCATE(BCKE_SUM) #endif END C --------------------- ASTASQ ------------------------------- SUBROUTINE ASTASQ(X,N,RI,ABA,KLP,SUS,SSQ,KLS) DIMENSION X(N,N) DOUBLE PRECISION ABA,SUS,SSQ C ESTIMATE AVERAGE OUTSIDE THE CIRCLE R = RI*RI NC = N/2+1 DO J=1,N T = J-NC XX = T*T DO I=1,N T=I-NC IF (XX + T * T .GT. R) THEN ABA = ABA + DBLE(X(I,J)) KLP = KLP + 1 ELSE SSQ = SSQ + X(I,J) * DBLE(X(I,J)) SUS = SUS + X(I,J) KLS = KLS + 1 ENDIF ENDDO ENDDO END C++********************************************************************* C--********************************************************************* SUBROUTINE FIXEDGE1(BCKP,NMAT,BCK3,N,IPCUBE,NN) DIMENSION BCKP(NMAT),BCK3(N,N,N),IPCUBE(5,NN) C------------------------------------------------------- C PUT ZEROS OUTSIDE C------------------------------------------------------- NT = 1 DO I=1,NMAT IF (NT .GT. NN) THEN BCKP(I) = 0.0 ELSEIF (I .LT. IPCUBE(1,NT)) THEN BCKP(I)= 0.0 ELSEIF(I .EQ. IPCUBE(2,NT)) THEN NT = NT+1 ENDIF ENDDO C ADD PIXELS ON THE EDGE C FIX THE EDGES IN BCKP DO K=1,N DO J=1,N DO I=1,N-1 IF (BCK3(I+1,J,K) .NE. 0.0) THEN BCK3(I,J,K) = BCK3(I+1,J,K) EXIT ENDIF ENDDO DO I=N,2,-1 IF (BCK3(I-1,J,K) .NE. 0.0) THEN BCK3(I,J,K) =BCK3(I-1,J,K) EXIT ENDIF ENDDO ENDDO ENDDO DO K=1,N DO I=1,N DO J=1,N-1 IF (BCK3(I,J+1,K) .NE. 0.0) THEN BCK3(I,J,K) = BCK3(I,J+1,K) EXIT ENDIF ENDDO DO J=N,2,-1 IF (BCK3(I,J-1,K) .NE. 0.0) THEN BCK3(I,J,K) = BCK3(I,J-1,K) EXIT ENDIF ENDDO ENDDO ENDDO DO J=1,N DO I=1,N DO K=1,N-1 IF (BCK3(I,J,K+1) .NE. 0.0) THEN BCK3(I,J,K) = BCK3(I,J,K+1) EXIT ENDIF ENDDO DO K=N,2,-1 IF (BCK3(I,J,K-1) .NE. 0.0) THEN BCK3(I,J,K) = BCK3(I,J,K-1) EXIT ENDIF ENDDO ENDDO ENDDO END C++********************************************************************* C--********************************************************************* SUBROUTINE FIXEDGE2(BCKP,NMAT,BCK3,N,IPCUBE,NN) DIMENSION BCKP(NMAT),BCK3(N,N,N),IPCUBE(5,NN) C------------------------------------------------------- C PUT ZEROS OUTSIDE C------------------------------------------------------- NT = 1 DO I=1,NMAT IF (NT .GT. NN) THEN BCKP(I) = 0.0 ELSEIF (I .LT. IPCUBE(1,NT)) THEN BCKP(I) = 0.0 ELSEIF (I .EQ. IPCUBE(2,NT)) THEN NT = NT+1 ENDIF ENDDO C ADD PIXELS ON THE EDGE C FIX THE EDGES IN BCKP DO K=1,N DO J=1,N DO I=2,N-1 IF (BCK3(I+1,J,K) .NE. 0.0) THEN BCK3(I,J,K) = BCK3(I+1,J,K) BCK3(I-1,J,K) = BCK3(I+1,J,K) EXIT ENDIF ENDDO DO I=N-1,2,-1 IF (BCK3(I-1,J,K) .NE. 0.0) THEN BCK3(I,J,K) = BCK3(I-1,J,K) BCK3(I+1,J,K) = BCK3(I-1,J,K) EXIT ENDIF ENDDO ENDDO ENDDO DO K=1,N DO I=1,N DO J=2,N-1 IF(BCK3(I,J+1,K).NE.0.0) THEN BCK3(I,J,K)=BCK3(I,J+1,K) BCK3(I,J-1,K)=BCK3(I,J+1,K) EXIT ENDIF ENDDO DO J=N-1,2,-1 IF(BCK3(I,J-1,K).NE.0.0) THEN BCK3(I,J,K) = BCK3(I,J-1,K) BCK3(I,J+1,K) = BCK3(I,J-1,K) EXIT ENDIF ENDDO ENDDO ENDDO DO J=1,N DO I=1,N DO K=2,N-1 IF(BCK3(I,J,K+1).NE.0.0) THEN BCK3(I,J,K) = BCK3(I,J,K+1) BCK3(I,J,K-1) = BCK3(I,J,K+1) EXIT ENDIF ENDDO DO K=N-1,2,-1 IF(BCK3(I,J,K-1).NE.0.0) THEN BCK3(I,J,K) = BCK3(I,J,K-1) BCK3(I,J,K+1) = BCK3(I,J,K-1) EXIT ENDIF ENDDO ENDDO ENDDO END C++********************************************************************* C--********************************************************************* SUBROUTINE FIXEDGE3(BCKP,NMAT,BCK3,N,IPCUBE,NN) DIMENSION BCKP(NMAT),BCK3(N,N,N),IPCUBE(5,NN) C------------------------------------------------------- C PUT ZEROS OUTSIDE C------------------------------------------------------- NT = 1 DO I=1,NMAT IF (NT.GT.NN) THEN BCKP(I) = 0.0 ELSEIF (I .LT. IPCUBE(1,NT)) THEN BCKP(I) = 0.0 ELSEIF (I .EQ. IPCUBE(2,NT)) THEN N T= NT+1 ENDIF ENDDO C ADD PIXELS ON THE EDGE C FIX THE EDGES IN BCKP DO K=1,N DO J=1,N DO I=3,N-1 IF(BCK3(I+1,J,K).NE.0.0) THEN BCK3(I,J,K)=BCK3(I+1,J,K) BCK3(I-1,J,K)=BCK3(I+1,J,K) BCK3(I-2,J,K)=BCK3(I+1,J,K) EXIT ENDIF ENDDO DO I=N-2,2,-1 IF(BCK3(I-1,J,K).NE.0.0) THEN BCK3(I,J,K)=BCK3(I-1,J,K) BCK3(I+1,J,K)=BCK3(I-1,J,K) BCK3(I+2,J,K)=BCK3(I-1,J,K) EXIT ENDIF ENDDO ENDDO ENDDO DO K=1,N DO I=1,N DO J=3,N-1 IF(BCK3(I,J+1,K).NE.0.0) THEN BCK3(I,J,K)=BCK3(I,J+1,K) BCK3(I,J-1,K)=BCK3(I,J+1,K) BCK3(I,J-2,K)=BCK3(I,J+1,K) EXIT ENDIF ENDDO DO J=N-2,2,-1 IF(BCK3(I,J-1,K).NE.0.0) THEN BCK3(I,J,K)=BCK3(I,J-1,K) BCK3(I,J+1,K)=BCK3(I,J-1,K) BCK3(I,J+2,K)=BCK3(I,J-1,K) EXIT ENDIF ENDDO ENDDO ENDDO DO J=1,N DO I=1,N DO K=3,N-1 IF(BCK3(I,J,K+1).NE.0.0) THEN BCK3(I,J,K)=BCK3(I,J,K+1) BCK3(I,J,K-1)=BCK3(I,J,K+1) BCK3(I,J,K-2)=BCK3(I,J,K+1) EXIT ENDIF ENDDO DO K=N-2,2,-1 IF(BCK3(I,J,K-1).NE.0.0) THEN BCK3(I,J,K)=BCK3(I,J,K-1) BCK3(I,J,K+1)=BCK3(I,J,K-1) BCK3(I,J,K+2)=BCK3(I,J,K-1) EXIT ENDIF ENDDO ENDDO ENDDO END C++********************************************************************* C--********************************************************************* SUBROUTINE BFIRSTS(BCKE,BCKP,N,ALA,IPCUBE,NN) DIMENSION BCKE(N,N,N),BCKP(N,N,N),IPCUBE(5,NN) c$omp parallel do private(kn,j,k,i) DO KN=1,NN J = IPCUBE(4,KN) K = IPCUBE(5,KN) DO I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN) BCKE(I,J,K) = BCKE(I,J,K)+ALA*(6*BCKP(I,J,K)-( & +BCKP(I+1,J,K)+BCKP(I,J+1,K)+BCKP(I,J,K+1) & +BCKP(I-1,J,K)+BCKP(I,J-1,K)+BCKP(I,J,K-1))) ENDDO ENDDO END C++********************************************************************* C--********************************************************************* SUBROUTINE BSECOND(BCKE,BCKP,N,ALA,IPCUBE,NN) DIMENSION BCKE(N,N,N),BCKP(N,N,N),ipcube(5,nn) c$omp parallel do private(kn,j,k,i) DO KN=1,NN J = IPCUBE(4,KN) K = IPCUBE(5,KN) DO I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN) BCKE(I,J,K) = BCKE(I,J,K) + ALA*(18*BCKP(I,J,K) & -4.0*BCKP(I+1,J,K)+BCKP(I+2,J,K) & -4.0*BCKP(I-1,J,K)+BCKP(I-2,J,K) & -4.0*BCKP(I,J+1,K)+BCKP(I,J+2,K) & -4.0*BCKP(I,J-1,K)+BCKP(I,J-2,K) & -4.0*BCKP(I,J,K+1)+BCKP(I,J,K+2) & -4.0*BCKP(I,J,K-1)+BCKP(I,J,K-2) & ) ENDDO ENDDO END C++********************************************************************* C--********************************************************************* SUBROUTINE BTHIRD(BCKE,BCKP,N,ALA,IPCUBE,NN) DIMENSION BCKE(N,N,N),BCKP(N,N,N),ipcube(5,nn) c$omp parallel do private(kn,j,k,i) DO KN=1,NN J = IPCUBE(4,KN) K = IPCUBE(5,KN) DO I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN) BCKE(I,J,K) = BCKE(I,J,K) + ALA*(60*BCKP(I,J,K) & -15.0*BCKP(I+1,J,K)+6.0*BCKP(I+2,J,K)-BCKP(I+3,J,K) & -15.0*BCKP(I-1,J,K)+6.0*BCKP(I-2,J,K)-BCKP(I-3,J,K) & -15.0*BCKP(I,J+1,K)+6.0*BCKP(I,J+2,K)-BCKP(I,J+3,K) & -15.0*BCKP(I,J-1,K)+6.0*BCKP(I,J-2,K)-BCKP(I,J-3,K) & -15.0*BCKP(I,J,K+1)+6.0*BCKP(I,J,K+2)-BCKP(I,J,K+3) & -15.0*BCKP(I,J,K-1)+6.0*BCKP(I,J,K-2)-BCKP(I,J,K-3) & ) ENDDO ENDDO END