! ! PSTG3R.SPLIT ! ! Version 1.007 30/06/21 NRB/CPB ! Version 0.999 01/07/13 CPB ! ! PSTG3R.SPLIT is a natural extension of the pstg3r code. ! If processors are available then diagonalise every Hamiltonian ! concurrently. Even if not, can carry-out multiple ones concurrently ! . ! ! Namelist : &prediag npw_per_subworld=1 icomplete_hdat=0 / ! ! npw_per_subworld=number of partial waves to diagonalize per subworld ! ! (npw/npw_per_subworld)*npcol*nprow MUST equal TOTAL number of procs, ! or simply, (number of sub worlds)*nprow*npcol. ! ! _____________________________________________________________________ ! ! Example : case with 6 partial waves, then use of ! ! &prediag npw_subworld=2 / ! (i.e. 2 partial waves per subworld) ! gives 3 subworlds, and use of ! ! &matrixdat npcol=2 nprow=2 / ! ! needs 3*2*2 processors, and so ! ! mpirun -np 12 ./pstg3r_split.x ! _____________________________________________________________________ ! ! ! icomplete_hdat = 1 will write header information to every H.DATXXX ! making it a 'complete' H.DAT in its own right. ! Useful for debugging or isolating a partial wave. ! ! icomplete_hdat = 0 default, then ! cat H.DAT[0-9][0-9][0-9] > H.DAT ! or use a merging code (e.g. hmerge) ! ! Do set IONEONE in NAMELIST STG3B to the number of symmetries in ! the STG2HXXX/RECUPHXXX files. (But unused if pstg2r_split files.) ! Then ! IONEONE=1 ! And ! npw_per_subworld = 1 ! is the fastest way to proceed as each subworld of processors ! just reads and diagonalizes the single STG2HXXX/RECUPHXXX ! Hamiltonian file it requires ! if ! IONEONE=0 (cluster default) then skips its way thru every file, ! IONEONE>0 then each subworld just goes to the file(s) it needs. ! ! Note: if ! npw_per_subworld = number of partial waves in sizeH.dat/sizeBP.dat ! then the run reduces to that of standard pstg3r ! ! ! ********************************************************************** ! 1.000 updates ! ! Cluster friendly (excludes pstg2r_split operation): ! ! IONEONE reactivated as the number of symmetries in STG2HXXX/RECUPHXXX ! and de-coupled from the number of subworlds (there was an ! implict one-to-one mapping of files and subworlds.) ! This enables subworlds to focus on only the files they need. ! ! zero-channel symmetries (inc. empty files) handled transparently ! (if npw_per_subworld=1 then those subworld processors do nothing, ! but it is typically one in a hundred.) ! ! ********************************************************************** ! 0.999 updates ! subroutine nochan : removes no channel partial waves ! from sizeH.dat files ! : (in case of npw_per_subworld=1 ) ! deletes empty STG2H files ! default: call nochan() commented out ! ! ********************************************************************** ! 0.99 updates ! ! a) Transparently uses pstg2r code that splits over target L S Pi. ! ! b) Provided DISTRIBUTE.DAT is present , you can excercise complete ! control over the number of processors given to each subworld ! ! based upon (pw_a)**3 / sum over (i (pw_i)**3 scaling) ! ! utility code : redistribute.f creates DISTRIBUTE.DAT ! !*********************************************************************** module comm_interface implicit none include "mpif.h" public comm_init ! Initialize MPI public comm_barrier ! MPI barrier public comm_finalize ! Terminate MPI public comm_split ! Split worlds integer, public :: icolour integer, public :: masteriam integer, public :: masternproc integer, public :: iam,nproc integer, public :: new_nodes integer, public :: new_comm integer, public :: npw,npw_per_subworld,icomplete_hdat integer, public :: nbp=0,nls=0 integer, public :: istart1 SAVE private integer :: mpicom CONTAINS !----------------------------------------------------------------------- subroutine comm_init() implicit none integer :: ier mpicom = MPI_COMM_WORLD call mpi_init(ier) call mpi_comm_rank(mpicom, masteriam, ier) call mpi_comm_size(mpicom, masternproc, ier) return end subroutine comm_init !----------------------------------------------------------------------- subroutine comm_barrier() implicit none integer :: ier call mpi_barrier(mpicom, ier) return end subroutine comm_barrier !----------------------------------------------------------------------- subroutine comm_finalize() implicit none integer :: ier call mpi_finalize(ier) return end subroutine comm_finalize !----------------------------------------------------------------------- subroutine comm_split(icolour) implicit none integer :: ier integer :: icolour call mpi_comm_split(mpicom,icolour,masteriam,new_comm,ier) call mpi_comm_rank(new_comm,iam,ier) call mpi_comm_size(new_comm,nproc,ier) return end subroutine comm_split !----------------------------------------------------------------------- end module comm_interface !*********************************************************************** program pstg3rsplit use comm_interface implicit real*8(a-h,o-z) logical ex,exdis include 'mpif.h' NAMELIST/PREDIAG/npw_per_subworld,icomplete_hdat C C C call comm_init() C C INQUIRE (FILE='dstg3',EXIST=EX) if(.NOT.EX)then continue else open(110,file='dstg3',status='unknown',form='formatted') icomplete_hdat=0 npw_per_subworld=1 read(110,prediag) close(110) endif C C call comm_barrier() C C DECIDE WHETHER WE ATE DOING A BP or LS calculation C C READ either sizeH.dat or sizeBP.dat to determine number of C partial waves C C INQUIRE (FILE='sizeBP.dat',EXIST=EX) if(.NOT.EX)then continue else open(111,file='sizeBP.dat',status='unknown',form='formatted') nbp=1 KOUNT=0 80 read(111,*,end=100)inumber KOUNT=KOUNT+1 goto 80 endif C INQUIRE (FILE='sizeH.dat',EXIST=EX) if(.NOT.EX)then continue else open(111,file='sizeH.dat',status='unknown',form='formatted') nls=1 KOUNT=0 90 read(111,92,end=100)char1 KOUNT=KOUNT+1 goto 90 endif C 92 FORMAT(A8) C if(iam.eq.0)write(0,*)'neither sizeH.dat or sizeBP.dat exist' C go to 1000 C 100 continue close(111) C npw=kount num_subworlds=npw/npw_per_subworld c if(num_subworlds*npw_per_subworld.ne.npw)then if(iam.eq.0)write(0,*)'npw/npw_per_subworld.ne.integer!' go to 1000 endif C C call comm_barrier() C C call nochan() C C call comm_barrier() C C icolour=mod(masteriam,npw) ! better to have procs sequential C INQUIRE(file='DISTRIBUTE.DAT',exist=exdis) if(.not.exdis)then num_proc_per_subworld=masternproc/num_subworlds if(num_proc_per_subworld*num_subworlds.ne.masternproc)then if(iam.eq.0) x write(0,*)'total number procs/num_subworlds.ne.integer!' go to 1000 endif icolour=masteriam/num_proc_per_subworld !standard else OPEN(444,file='DISTRIBUTE.DAT',form='formatted',status='unknown') DO I=0,NPW-1 READ(444,*)istart1,ifinish1 if((masteriam.ge.istart1).and.(masteriam.le.ifinish1))icolour=I ENDDO CLOSE(444) endif call comm_barrier() !don't remove call comm_split(icolour) C C C izero=-1 if(iam.eq.0)izero=icolour call MPI_Bcast(izero,1,mpi_integer,0,new_comm,ier) C write(0,*)'old id= ',masteriam,'new id =',iam,'icolour= ',icolour C call pstg3rsub() C call comm_barrier() C if(masteriam.eq.0)then call cpu_time(timef) time=timef-timei C write(0,999) time/60.,nproc 999 format(//1x,'CPU TIME=',f9.3,' MIN -- processors=:',i4) C endif C 1000 continue call comm_barrier() call comm_finalize() stop C C end C*********************************************************************** c PSTG3RSUB :: derived from c PSTG3R :: Parallel Version 4.4 of stg3r. 08/08/13 c c c D. M. Mitnik c Version 1.0: c Subroutines DMAT, DMAT1 and DMAT2 deleted c Allows use in Breit-Pauli mode c Version 2.0: every processor reads , without broadcasting, c except for UNIT5. c Version 3.2: test dimensions option c Version 3.4: Interface with parallel stg2r STG2HXXX.DAT c Version 3.8: Rework Allocate c c C. P. Ballance c Version 3.9: Rework gathering of surface amplitudes. c Interface with Parallel pstgjk/precupd. c INTEGER*8 for large cases. c Version 4.1: Partitioned R-matrix. c Version 4.3: block writing of wmat (NDIV in matrixdat) c c N. R. Badnell c Version 2.1: Including an alternative use of pdsyevd c Version 3.1: Memory allocatable. c Version 3.3: Allocation mods to HAMNEWx suite c Version 3.5: Use single UNIT number for STG2HXXX.DAT files c Version 3.6: Remove more broadcasting, let each proc compute. c Version 3.7: Multiple AMATUXXX.DAT for parallel pstg2r. c Version 4.1: Partitioned R-matrix. c Version 4.2: KAB's method for eliminating weakly coupled N+1 terms c Version 4.3: block writing of wmat (NDIV in matrixdat) c Version 4.4: IONEONE skips file reads for IPWINT>1. c C*********************************************************************** C C N. R. BADNELL UoS v3.4 30/05/08 C C C*********************************************************************** C C Belfast Originated Atomic R-matrix Codes C C*********************************************************************** C C THE THIRD PART OF C C A GENERAL PROGRAM TO CALCULATE ATOMIC CONTINUUM C C PROCESSES USING THE R-MATRIX METHOD C C S T G 3 C C ORIGINALLY DISTRIBUTED BY C C QUEEN'S UNIVERSITY BELFAST C C*********************************************************************** C C DIAGONALIZES THE HAMILTONIANS AND PROCESSES DIPOLE MATRIX C ELEMENTS USING AS INPUT THE H (AND PERHAPS THE D) FILES C CREATED IN STG2 OR RECUPD - dipole operation not coded. C C CAN ALSO HANDLE H AND D FILES PRODUCED BY DARC, & DTO3 INTERFACE. C C WILL USE A PARTIONED R-MATRIX IF IPRCENT SPECIFIED (40-100). C C USES SCALAPACK PACKAGE FOR DIAGONALIZING. C C*********************************************************************** C C ROUTINES USED IN STG3 C C*********************************************************************** C C MNSTG3 C STG3 DRIVER C DA2 C DIAG C DMAT !! parallel !! DMAT not implemented C DMAT1 C DMAT2 C HAMNEW0 C HAMNEW C ORDER C RECOV2 C RSCT C STG3RD C TAPERD C VMUL C WRITOP c..................**** parallel ****.......................... c initctxt ! initializes the BLACS c read2r ! reads in "sizeH.dat" the size of H c Hsize ! calculates the size of the local H c locH ! mapping the global H into the local H c diag_pdsyevd ! driver for pdsyevd SCAlapack diagonalization c diag_pdsyevx ! driver for pdsyevx SCAlapack diagonalization c demapping ! demapping the local position (i,j) into (I,J) c..................**** parallel ****.......................... C C*********************************************************************** C C INPUT/OUTPUT CHANNELS USED IN STG3 C C*********************************************************************** C C IREAD USER INPUT FILE (FILE dstg3) C IWRITE OUTPUT TO LINE PRINTER (FILE rout3r) C C IDISC1 OPTIONAL SCRATCH FILE FOR THE H-MATRIX C IDISC2 SCRATCH DIRECT ACCESS FILE OF H EIGENVECTORS C IDISC3 SCRATCH DIRECT ACCESS FILE OF PARTIALLY C PROCESSED DIPOLE MATRIX ELEMENTS C C ITAPE1 INPUT FILE OF STG2 DIPOLE MATRICES C ITAPE2 INPUT FILE OF STG2 H-MATRICES C ITAPE3 OUTPUT FILE OF DIAGONALIZED H-MATRICES C ITAPE4 OUTPUT FILE OF PROCESSED DIPOLE MATRICES. C c **** parallel **** C IREAD (5) .. parallel input data .. dstg3 c **** parallel **** C IWRITE (6) .. printed output .. rout3r C C IPUNCH .. NOT USED C C IDISC1 (11) .. scratch C IDISC2 (12) .. scratch (DA2) C IDISC3 (13) .. scratch C IDISC4 .. NOT USED C C ITAPE1 (1) .. RECUPD/STG2 dump (dipole) .. RECUPD/STG2D.DAT C if IPOLPH=2 C ITAPE2 (2) .. RECUPD/STG2 dump (hamiltonians) .. RECUPH/STG2H.DAT C always used C ITAPE3 (3) .. STG3 dump (hamiltonians) .. H.DAT C always used C ITAPE4 (4) .. STG3 dump (dipole matrices) .. D00 etc C if IPOLPH=2 C c **** parallel **** c unit(30) .. size of Hamiltonians .. sizeH.dat (input) c unit(35) .. CPU time information .. time**.dat (output) c **** parallel **** C C*********************************************************************** C C DIMENSIONING PARAMETERS USED IN PSTG3 C C*********************************************************************** C C INCLUDE PARAMETERS: C C CHF (75) HIGHEST NUMBER OF CHANNELS (NCHAN) C LMX (8) MULTIPOLES IN POTENTIAL (LAMAX) C LR2 (20) HIGHEST L+1 FOR CONTINUUM ORBITALS (LRANG2) C MEG (1) MEGA-WORDS OF MEMORY TO REDUCE I/O (CAN BE 0) CFIXED KIL (100) KILO-WORDS OF MEMORY TO REDUCE I/O (CAN BE 0) C NC2 (600) N+1 ELECTRON CONFIGURATIONS FOR GIVEN S L PI (NCFGP) C NR2 (40) NUMBER OF CONTINUUM ORBITALS FOR GIVEN L (NRANG2) C SLP (80) NUMBER OF DIFFERENT N+1 ELECTRON SYMMETRIES (INAST) C TAR (63) TARGET STATES TOTAL (NAST) C C*********************************************************************** C C MODULES USED IN PSTG3 C C*********************************************************************** C c param c basic2 c boundH c HAMNUW C C*********************************************************************** module param implicit real*8(a-h,o-z) include 'PARAM' end c*********************************************************************** module basic2 use param implicit real*8(a-h,o-z) allocatable cf(:,:,:) dimension et(mzchf),kj(mzchf),l2p(mzchf),lstarg(mzchf) x ,nconat(mztar) integer lrgl,nspn,npty end c*********************************************************************** module boundH use param implicit real*8(a-h,o-z) allocatable hbc(:,:),hbb(:,:) end c*********************************************************************** module HAMNUW use param implicit real*8(a-h,o-z) allocatable a(:,:),xo(:,:),xl(:),work(:),b(:,:) end c*********************************************************************** c subroutine pstg3rsub() use param use comm_interface IMPLICIT REAL*8 (A-H,O-Z) C C PARAMETER (MXMEM=MZMEG*1000000+MZKIL*1000+1) PARAMETER (MXN21=MZNR2+1) PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) c PARAMETER (MXTRI=(MXMAT*(MXMAT+1))/2) C C MXDM=MAX(MZNC2,MZNR2): C c PARAMETER (MXD1=MZNC2/MZNR2,MXD2=MZNR2/MZNC2,MXD3=MXD1+MXD2, c A MXDM=MZNC2*MXD1/MXD3+MZNR2*MXD2/MXD3) C C MXMATX=MAX(MXDMAT,MXRSCT)+1 = SIZE OF /MATRIX/ IN DMAT1, RSCT: C c PARAMETER (MXDMAT=4*MXMAT*MZCHF+2*MXMAT*MXDM+2*MZNR2*MXDM, c A MXRSCT=MXTRI+MZCHF*MXMAT+MZNR2*MXMAT) c PARAMETER (MXM1=MXDMAT/MXRSCT,MXM2=MXRSCT/MXDMAT,MXM3=MXM1+MXM2, c A MXMATX=MXDMAT*MXM1/MXM3+MXRSCT*MXM2/MXM3+1) C C*********************************************************************** C C COMMON BLOCKS USED IN STG3 C C*********************************************************************** C COMMON /ADDE/EST(MZTAR),EPOLE,IPOLE,NEST,AST(MZTAR),MERGE(MZTAR) COMMON /ALPHA/LSP(MZSLP,3),LCH(MZSLP),LCFG(MZSLP) COMMON /BASIC1/BSTO,RA,NELC,NZ,LRANG2,NRANG2,MAXNHF(MZLR2), A MAXNLG(MZLR2),MAXNC(MZLR2),LRANG1 COMMON /BASIN/EIGENS(MZNR2,MZLR2),ENDS(MXN21,MZLR2),DELTA,ETA COMMON /BUTT/COEFF(3,MZLR2),EK2MAX,EK2MIN,MAXNCB(MZLR2),NELCOR COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST,IONEONE COMMON /COPY/ICOPY1,ICOPY2,ITOTAL,ICOUNT COMMON /CUPINT/MNP1,NCONHP,NCHAN,N2HDAT,NTARSYM COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /INFORM/IREAD,IWRITE,IPUNCH c COMMON /MATRIX/DUME(MXMATX) !! parallel !! DMAT deleted c COMMON /MEMORY/ARRAY(MXMEM),MEM1,MREC1 COMMON /MORDER/LAMAX,LAMBC,LAMCC,LAMIND,LAM COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /NRBDIP/MAXLD COMMON /RCASES/NOT1,NOT2,NOT3,NPROG,NCFGP COMMON /RECORD/ILSPI(4,MZSLP),IORDER(MZCHF,MZSLP) COMMON /RECOV/IPLACE COMMON /REDMEL/CGC(20),MAXM1 C COMMON /SPACES/XX(MZNR2,MXMAT) COMMON /STATE/ENAT(MZTAR),LAT(MZTAR),ISAT(MZTAR),LPTY(MZTAR),NAST C C c **** parallel **** include 'mpif.h' common /bksize/nl,nr common /contsize/icontold,nonzer common /crashblock/istop common /globsize/idimH(mzslp) C common /parablock/iam,nproc common /procdat/ictxt,myrow,mycol common /matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest X ,NDIV c c call mpi_init(ierr) c call mpi_comm_rank(mpi_comm_world,iam,ierr) c call mpi_comm_size(mpi_comm_world,nproc,ierr) istop=0 c **** parallel **** C C C*********************************************************************** C C C STG3 MAIN PROGRAM C C*********************************************************************** C C MEM1 AND MREC1 ARE THE MEMORY AND DA FILE POINTERS C CMICD GETCPUS MEM1 = 0 MREC1 = -1 call cpu_time(timei) c write(0,*) ' begin proc=:',iam c c c CALL STG3 C C if (nproctest.ne.0) go to 1000 call blacs_gridexit(ictxt) call comm_barrier() C call mpi_comm_free(new_comm,ier) C if (istop.ne.0) then call mpi_barrier(mpi_comm_world,ierr) go to 1000 endif C c **** parallel **** if (iam.eq.0) then call cpu_time(timef) time=timef-timei time=time/60. write(35,2999) time 2999 format(//2x,' Total CPU time=:',1pg10.3,' min.') write(iwrite,999) time,nproc 999 format(//1x,'CPU TIME=',f9.3,' MIN -- processors=:',i4) endif 1000 continue close(35) c **** parallel **** c call comm_barrier() c call comm_finalize() c STOP END C*********************************************************************** SUBROUTINE DIAG(N,IUP,Z,D,E,MXMAT) C use param IMPLICIT REAL*8 (A-H,O-Z) C C BADNELL & BURGESS D.A.M.T.P. CAMBRIDGE C C DIAGONALIZATION OF REAL SYMMETRIC N-BY-N MATRIX Z. C C METHOD: HOUSEHOLDER REDUCTION TO TRI-DIAGONAL FORM AND SHIFTED C QL ALGORITHM TO DETERMINE THE E-VALUES AND E-VECTORS. C C BASED ON MARTIN, REINSCH & WILKINSON: NUM. MATH. 11, 181-95 (1968). C C INPUT REQUIRED. N, IUP AND Z. ONLY LOWER TRIANGLE OF Z NEED BE SUPPLIED. C MATRIX Z OVERWRITTEN BY EIGENVECTORS OF Z. C IUP=1/-1 ASC/DESCENDING SORT, 0 NO SORT. C MXMAT, IS THE ROW DIMENSION OF Z IN THE CALLING ROUTINE. C C OUTPUT. Z AND D, WHERE Z CONSISTS OF COLUMN EIGENVECTORS C AND D CONSISTS OF CORRESPONDING EIGENVALUES. C C NOTE: E IS A WORKING ARRAY. C PARAMETER (TOL = 1.0D-75) PARAMETER (EPS = 1.0D-15) PARAMETER (ZERO = 0.0D0) PARAMETER (ONE = 1.0D0) PARAMETER (JMAX = 30) C DIMENSION D(N),E(N),Z(MXMAT,N) C COMMON /INFORM/IREAD,IWRITE,IPUNCH C DO 1 I = 1,N D(I) = Z(N,I) 1 CONTINUE IF (N.LE.1) GO TO 20 C C HOUSEHOLDER REDUCTION TO TRI-DIAGONAL FORM C DO 19 I = N,2,-1 L = I - 1 F = D(I-1) G = ZERO DO 5 K = 1,I-2 G = G + D(K)*D(K) 5 CONTINUE H = G + F*F IF (G.GT.TOL) GO TO 8 E(I) = F H = ZERO DO 7 J = 1,L D(J) = Z(L,J) Z(I,J) = ZERO Z(J,I) = ZERO 7 CONTINUE GO TO 18 8 G = SQRT(H) IF (F.GE.ZERO) G = -G E(I) = G H = H - F*G D(L) = F - G DO 14 J = 1,L E(J) = ZERO 14 CONTINUE DO 15 J = 1,L Z(J,I) = D(J) G = E(J) + Z(J,J)*D(J) DO 13 K = J+1,L G = G + Z(K,J)*D(K) E(K) = E(K) + Z(K,J)*D(J) 13 CONTINUE E(J) = G 15 CONTINUE F = ZERO DO 12 J = 1,L E(J) = E(J)/H F = F + E(J)*D(J) 12 CONTINUE HH = F/(H+H) DO 11 J = 1,L E(J) = E(J) - HH*D(J) 11 CONTINUE DO 17 J = 1,L F = D(J) G = E(J) DO 16 K = J,L Z(K,J) = Z(K,J) - F*E(K) - G*D(K) 16 CONTINUE D(J) = Z(L,J) Z(I,J) = ZERO 17 CONTINUE 18 D(I) = H 19 CONTINUE C C C ACCUMULATE TRANSFORMATION MATRICES C DO 28 I = 2,N L = I - 1 Z(N,L) = Z(L,L) Z(L,L) = ONE H = D(I) IF (H.EQ.ZERO) GO TO 25 DO 21 K = 1,L D(K) = Z(K,I)/H 21 CONTINUE DO 24 J = 1,L G = ZERO DO 22 K = 1,L G = G + Z(K,I)*Z(K,J) 22 CONTINUE DO 23 K = 1,L Z(K,J) = Z(K,J) - G*D(K) 23 CONTINUE 24 CONTINUE 25 DO 27 J = 1,L Z(J,I) = ZERO 27 CONTINUE 28 CONTINUE DO 29 I = 1,N D(I) = Z(N,I) Z(N,I) = ZERO 29 CONTINUE 20 E(1) = ZERO Z(N,N) = ONE C C C SHIFTED QL ALGORITHM TO DETERMINE E-VALUES & E-VECTORS C DO 32 I = 2,N E(I-1) = E(I) 32 CONTINUE E(N) = ZERO B = ZERO F = ZERO DO 54 L = 1,N J = 0 H = EPS*(ABS(D(L))+ABS(E(L))) IF (B.LT.H) B = H DO 36 M = L,N IF (ABS(E(M)).LE.B) GO TO 37 36 CONTINUE 37 IF (M.EQ.L) GO TO 53 38 IF (J.EQ.JMAX) GO TO 62 J = J+ 1 P = E(L) + E(L) G = D(L) H = D(L+1) - G IF (ABS(H).GE.ABS(E(L))) GO TO 43 P = H/P R = SQRT(P*P+ONE) H = P + R IF (P.LT.ZERO) H = P - R D(L) = E(L)/H GO TO 44 43 P = P/H R = SQRT(P*P+ONE) D(L) = E(L)*P/(R+ONE) 44 H = G - D(L) DO 46 I = L+1,N D(I) = D(I) - H 46 CONTINUE F = F + H P = D(M) C = ONE S = ZERO DO 52 I = M-1,L,-1 G = C*E(I) H = C*P IF (ABS(P).LT.ABS(E(I))) GO TO 49 C = E(I)/P R = SQRT(C*C+ONE) E(I+1) = S*P*R S = C/R C = ONE/R GO TO 50 49 C = P/E(I) R = SQRT(C*C+ONE) E(I+1) = S*E(I)*R S = ONE/R C = C/R 50 P = C*D(I) - S*G D(I+1) = H + S*(C*G+S*D(I)) DO 51 K = 1,N H = Z(K,I+1) Z(K,I+1) = S*Z(K,I) + C*H Z(K,I) = C*Z(K,I) - S*H 51 CONTINUE 52 CONTINUE E(L) = S*P D(L) = C*P IF (ABS(E(L)).GT.B) GO TO 38 53 D(L) = D(L) + F 54 CONTINUE C IF(IUP.EQ.0)RETURN C C BEGIN SORTING INTO ASCENDING E-VALUES C DO 61 I = 1,N K = I P = D(I) DO 57 J = I+1,N IF (IUP.GT.0.AND.D(J).GT.P) GO TO 57 IF (IUP.LT.0.AND.D(J).LT.P) GO TO 57 K = J P = D(J) 57 CONTINUE IF (K.EQ.I) GO TO 61 D(K) = D(I) D(I) = P DO 60 J = 1,N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 60 CONTINUE 61 CONTINUE C RETURN C C 62 WRITE(IWRITE,100) 100 FORMAT(' FAILED IN DIAG, TOO MANY ITERATIONS') RETURN C END C********************************************************************** SUBROUTINE HAMNEW0 C use param use basic2 use HAMNUW use comm_interface IMPLICIT REAL*8 (A-H,O-Z) C C----------------------------------------------------------------------- C c **** dario **** C THIS ROUTINE READS IN THE OVERLAP MATRIX FROM STG2 c **** dario **** C C----------------------------------------------------------------------- C LOGICAL EXA CHARACTER*12 IFLNM C COMMON /REDTOL/TOLER COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST,IONEONE COMMON /INFORM/IREAD,IWRITE,IPUNCH C SAVE IFL DATA IFL/1/ c **** parallel **** include 'mpif.h' common /crashblock/istop C common /parablock/iam,nproc common /globsize/idimH(mzslp) common /bksize/nl,nr common /contsize/icontold,nonzer character*1 filec,filed,fileu c **** parallel **** C C C C RETURN (DO NOTHING TO HAMILTONIAN) IF TOLERANCE IS REAL SMALL C IF(TOLER.LT.1.D-19)RETURN C if(inast.ne.0)rewind(77) C c1 READ(77,*,END=3)NSPIN,NLAM,NPAR,NL,NR,NONZER 1 READ(77,END=3)NSPIN,NLAM,NPAR,NL,NR,NONZER c IF((NSPIN.NE.NSPN.OR.NLAM.NE.LRGL.OR.NPAR.NE.NPTY) X .and.inast.ne.0)THEN !Parallel: SLP not read yet C DO N=1,NONZER c READ(77,*)I,J,DUM READ(77)I,J,DUM ENDDO GO TO 1 ENDIF C IF(NONZER.EQ.0)RETURN nt=max(nl,nr) allocate(A(nt,nt),XO(nt,nt),XL(nt),WORK(nt),B(nt,nt),stat=ierr) if(ierr.ne.0)stop 'allocation fails in HAMNEW0' C C INITIALIZE A C DO J=1,NR DO I=1,NL A(I,J)=0.0D0 ENDDO ENDDO icontold = idimH(mskip)-nr ! this is ncfgp2 C C READ IN NON-ZERO ELEMENTS OF A-MATRIX FROM STG2 C DO 10 NZ=1,NONZER c READ(77,*)I,J,ATEMP READ(77)I,J,ATEMP A(I,J)=A(I,J)+ATEMP 10 CONTINUE c if (mskip.lt.ipwinit) then deallocate(A,XO,XL,WORK,B) return endif C C MAKE XO=A*AT C DO I=1,NL DO J=1,NR B(J,I)=A(I,J) ENDDO ENDDO DO J=1,NL DO I=1,NL XO(I,J)=0.0D0 ENDDO DO K=1,NR DO I=1,NL XO(I,J)=XO(I,J)+A(I,K)*B(K,J) ENDDO ENDDO ENDDO C C MAKE O AND L=OT*A*AT*O (DIAGONAL EIGENVALUES) C CALL DIAG(NL,-1,XO,XL,WORK,nt) C if(iam.eq.0)then IF(NBUG9.GT.0)THEN c WRITE(78,*)' EIGENVALUES FOR PARTIAL WAVE ',NSPN2,LRGL2,NPTY2 WRITE(78,*)' EIGENVALUES FOR PARTIAL WAVE ',nspin,nlam,npar DO I=1,NL WRITE(78,100)XL(I) 100 FORMAT(6E13.6) ENDDO ENDIF endif C DO I=1,NL XL(I)=ABS(XL(I)) IF(XL(I).LT.1.D-18)XL(I)=1.D-18 ENDDO C C MAKE B=OT*A C DO J=1,NR DO I=1,NL SUM=0.0D0 DO K=1,NL SUM=SUM+XO(K,I)*A(K,J) ENDDO B(I,J)=SUM/SQRT(XL(I)) ENDDO ENDDO C C MAKE ANEW=(1/SQRT(L))*OT*A C DO J=1,NR DO I=1,NL A(I,J)=B(I,J) ENDDO ENDDO C C DETERMINE HOW MANY BASIS FUNCTIONS WILL BE KEPT C IF (NL.GT.NR)NL=NR NLNEW=NL DO I=1,NL IF(ABS(XL(NL-I+1)).LT.TOLER)THEN NLNEW=NLNEW-1 ELSE GO TO 51 ENDIF ENDDO c.....recalculate the Hamiltonian size 51 idiff=nr-nlnew ! channels omitted in new H idimH(mskip)=idimH(mskip)-idiff NL=NLNEW C C WRITE NEW MATRIX C if(iam.eq.0)then IF(NBUG9.GT.10)THEN WRITE(78,*)' NEW A-MATRIX OF DIMENSION ',NL,'X',NR DO J=1,NR WRITE(78,100)(A(I,J),I=1,NL) ENDDO ENDIF endif C return C C OPEN NEW FILE AMATUXXX.DAT 3 ich0 = ichar('0') ic = ifl/100 id = (ifl - 100*ic)/10 iu = (ifl - 100*ic - 10*id) ic = ic + ich0 id = id + ich0 iu = iu + ich0 filec = char(ic) filed = char(id) fileu = char(iu) IFLNM='AMATU'//filec//filed//fileu//'.DAT' INQUIRE (FILE=IFLNM,EXIST=EXA) IF(EXA)THEN CLOSE(77) OPEN (UNIT=77,FILE=IFLNM,STATUS='OLD',FORM='UNFORMATTED') WRITE(IWRITE,3210) IFLNM 3210 FORMAT(' OPENING FILE ',A12/) IFL=IFL+1 GO TO 1 ELSE WRITE(IWRITE,*)' SR.HAMNEW0: ERROR, NS,NL,NP =', C nspin,nlam,npar,' CASE NOT FOUND' if(iam.ne.0)write(0,*)' SR.HAMNEW0: ERROR, SLP CASE NOT FOUND' istop=1 return ENDIF 5555 format(10(1x,f9.4)) end C********************************************************************** SUBROUTINE HAMNEW(HNP1,irowsize,icolsize, + MNP1,NCFGP2,LRGL2,NSPN2,NPTY2) C use param use boundH use HAMNUW use comm_interface IMPLICIT REAL*8 (A-H,O-Z) C C----------------------------------------------------------------------- C C THIS ROUTINE TRANSFORMS THE HAMILTONIAN, PERHAPS REDUCING THE NUMBER OF C (N+1)-ELECTRON CONFIGURATIONS - TWG & NRB C C----------------------------------------------------------------------- C COMMON /REDTOL/TOLER COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /INFORM/IREAD,IWRITE,IPUNCH c **** parallel **** include 'mpif.h' common /crashblock/istop C common /parablock/iam,nproc common /globsize/idimH(mzslp) common /procdat/ictxt,myrow,mycol common /contsize/icontold,nonzer common /bksize/nl,nr dimension HNP1(irowsize,icolsize) c **** parallel **** C C RETURN (DO NOTHING TO HAMILTONIAN) IF TOLERANCE IS REAL SMALL C IF(TOLER.LT.1.D-19)RETURN C TRANSFORM (AND REDUCE) HAMILTONIAN C C HCC -> HCC (do nothing) C HBC -> A*HBC C HCB -> HCB*AT C HBB -> A*HBB*AT C NCONT=MNP1-NCFGP2 MNPN=NCONT+NL c DO J=1,NL !transpose A DO I=1,NCFGP2 XO(I,J)=A(J,I) ENDDO ENDDO c c xo=transpose(a) !could try f90 intrinsic...NRB C C C B-C C do 100 i=1,nl do 90 j=1,ncont c IH and JH are the I,J coordinates of H IH = i+ncont JH = j call mapping(JH,IH,irow,icol,ipr,ipc) if (ipr.ne.myrow.or.ipc.ne.mycol) go to 90 sumij = 0.0D0 do 50 k=1,ncfgp2 sumij = sumij + XO(k,i)*HBC(k,j) 50 continue HNP1(irow,icol) = sumij 90 continue 100 continue C C B-B C do j=1,nl do i=1,ncfgp2 B(i,j) = 0.0D0 enddo do k=1,ncfgp2 do i=1,ncfgp2 B(i,j) = B(i,j) + HBB(i,k)*XO(k,j) enddo enddo enddo do 300 i=1,nl do 290 j=i,nl c IH and JH are the I,J coordinates of H IH = i+ncont JH = j+ncont call mapping(IH,JH,irow,icol,ipr,ipc) if (ipr.ne.myrow.or.ipc.ne.mycol) go to 290 sumij = 0.0D0 do 250 k=1,ncfgp2 sumij = sumij + XO(k,i)*B(k,j) 250 continue HNP1(irow,icol) = sumij 290 continue 300 continue C C RE-DEFINE MNP1, NCFGP2 C IF(NL.LT.NCFGP2)THEN WRITE(IWRITE,*)' REDUCING # OF (N+1)-CONFIGS FROM ', + NCFGP2,' TO ',NL NCFGP2=NL MNP1=NCONT+NL ENDIF c deallocate(a,xo,xl,work,b) C RETURN C END C********************************************************************** SUBROUTINE HAMNEW2(HNP1,irowsize,icolsize, + MNP1,NCFGP2,LRGL2,NSPN2,NPTY2) C use param use boundH use HAMNUW use comm_interface IMPLICIT REAL*8 (A-H,O-Z) C C----------------------------------------------------------------------- C C THIS ROUTINE SHIFTS THE EIGENVALUES OF THE BOUND-BOUND C HAMILTONIAN - TWG C C----------------------------------------------------------------------- C COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /STATE/ENAT(MZTAR),LAT(MZTAR),ISAT(MZTAR),LPTY(MZTAR),NAST COMMON /SHIFT/ESHIFT,NSHIFT COMMON /INFORM/IREAD,IWRITE,IPUNCH C c **** parallel **** include 'mpif.h' common /crashblock/istop C common /parablock/iam,nproc common /globsize/idimH(mzslp) common /procdat/ictxt,myrow,mycol dimension HNP1(irowsize,icolsize) common /contsize/icontold,nonzer common /bksize/nl,nr c **** parallel **** C IF(NSHIFT.LE.0)RETURN C READ(81,*)NSPIN,NLAM,NPAR IF(NSPIN.NE.NSPN2.OR.NLAM.NE.LRGL2.OR.NPAR.NE.NPTY2)THEN WRITE(IWRITE,*)' ERROR, DIDNT FIND CORRECT PARTIAL WAVE INFO' WRITE(IWRITE,*)NSPIN,NLAM,NPAR,' ON FILE ',NSPN2,LRGL2,NPTY2, + ' PROGRAM' istop=1 return ENDIF c allocate(A(ncfgp2,ncfgp2),XL(ncfgp2) x,WORK(ncfgp2),stat=ierr) if(ierr.ne.0)stop'allocation fails HAMNEW2' C C C DIAGONALIZE THIS C CALL DIAG(NCFGP2,-1,hbb,XL,WORK,ncfgp2) C if (iam.eq.0)WRITE(79,*) C ' EIGENVALUES FOR PARTIAL WAVE ',NSPN2,LRGL2, CNPTY2,'RELATIVE TO GROUND CONTINUUM (IN RYD)' NBELOW=0 DO I=1,NCFGP2 EREL=2.0D0*(XL(I)-ENAT(1)) IF(EREL.LT.0.0D0)NBELOW=NBELOW+1 if (iam.eq.0)WRITE(79,111)I,EREL 111 FORMAT(I5,E13.6) ENDDO C C CORRECT EIGENVALUES C READ(81,*)NUM IF(NUM.LE.0) GO TO 81 DO I=1,NUM READ(81,*)NCORRECT,DELTA J=NCFGP2-(NCORRECT+NBELOW)+1 XL(J)=XL(J)+0.5D0*DELTA ENDDO C C MAKE NEW BOUND-BOUND HAMILTONIAN C do j=1,ncfgp2 do i=1,ncfgp2 a(i,j)=0.0d0 enddo enddo DO K=1,NCFGP2 DO J=1,NCFGP2 T=XL(K)*hbb(J,K) DO I=1,NCFGP2 A(I,J)=A(I,J)+hbb(I,K)*T ENDDO ENDDO ENDDO C C PUT HAMILTONIAN BACK ONTO HNP1 VECTOR C NCONT=MNP1-NCFGP2 C K=0 DO 80 I=1,NCFGP2 DO 82 J=I,NCFGP2 c IH and JH are the I,J coordinates of H IH = i+ncont JH = j+ncont call mapping(IH,JH,irow,icol,ipr,ipc) if (ipr.ne.myrow.or.ipc.ne.mycol) go to 82 HNP1(irow,icol) = A(I,J) 82 CONTINUE 80 CONTINUE C 81 deallocate(a,xl,work) c RETURN END C********************************************************************** SUBROUTINE ORDER(EN,NORDER,NDIM,IUP) C use param IMPLICIT REAL*8 (A-H,O-Z) C C C C----------------------------------------------------------------------- C C RETURNS NORDER(I)=POINTER TO I-TH ENERGY IN EN ARRAY, C IUP=1 FOR ASCENDING ENERGIES, IUP=-1 FOR DESCENDING ENERGIES C C----------------------------------------------------------------------- DIMENSION EN(NDIM),NORDER(NDIM) C----------------------------------------------------------------------- DO 40 K = 1,NDIM J = K J1 = J - 1 IF (J1.EQ.0) GOTO 30 E = EN(J) + IUP*1.0D-7 DO 20 I = 1,J1 IF (J.LT.K) GOTO 10 N = NORDER(I) IF (IUP.GT.0 .AND. E.GT.EN(N)) GOTO 20 IF (IUP.LT.0 .AND. E.LT.EN(N)) GOTO 20 10 CONTINUE NORDER(J) = NORDER(J-1) J = J - 1 20 CONTINUE 30 NORDER(J) = K 40 CONTINUE C END C********************************************************************** SUBROUTINE RECOV2(SUBNAM,PARNAM,IDIM,NEED) C use param use comm_interface IMPLICIT REAL*8 (A-H,O-Z) C C C----------------------------------------------------------------------- C C THIS ROUTINE IS CALLED ONLY FOR ARRAY OVERFLOW C C SUBNAM = SUBROUTINE NAME C PARNAM = PARAMETER NAME C IDIM = CURRENT DIMENSION C NEED = REQUIRED DIMENSION, RETURN NEED=IDIM C IPLACE = 0 TO STOP PROGRAM, C OTHERWISE RETURN IPLACE=NEED TO THE CALLING ROUTINE. C C----------------------------------------------------------------------- CHARACTER*6 SUBNAM,PARNAM,PREPRO C COMMON /INFORM/IREAD,IWRITE,IPUNCH COMMON /RECOV/IPLACE c **** parallel **** include 'mpif.h' common /crashblock/istop c common /parablock/iam,nproc c **** parallel **** C----------------------------------------------------------------------- PREPRO = PARNAM IF (PREPRO(1:3).EQ.' MZ') PREPRO(1:3) = ' MZ' WRITE (IWRITE,3000) SUBNAM,PREPRO,IDIM,PREPRO,NEED C IF (IPLACE.EQ.0) THEN WRITE (IWRITE,3010) istop = 1 return ENDIF C IF (IPLACE.LT.0) WRITE (IWRITE,3020) IPLACE = NEED NEED = IDIM C/ 3000 FORMAT (/' * ARRAY OVERFLOW IN ', A A6/' MUST INCREASE DIMENSION GIVEN BY ',A6,' =',I7, B ' TO AT LEAST ',A6,' =',I9) 3010 FORMAT (/' PROGRAM TERMINATES IN RECOV2'/) 3020 FORMAT (/' CHECK TO SEE IF OTHER ARRAYS ARE GOING TO BE EXCEEDED') END C********************************************************************** SUBROUTINE STG3 c use param use basic2 use comm_interface IMPLICIT REAL*8 (A-H,O-Z) C C C C----------------------------------------------------------------------- LOGICAL NBUG10 C COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST,IONEONE COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /INFORM/IREAD,IWRITE,IPUNCH COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /RECOV/IPLACE COMMON /CUPINT/MNP1,NCONHP,NCHAN,N2HDAT,NTARSYM C LOGICAL EX,ESPLIT c **** parallel **** COMMON /REDTOL/TOLER COMMON /SHIFT/ESHIFT,NSHIFT common /crashblock/istop include 'mpif.h' allocatable :: HNP1(:,:) character*1 :: filecc,filedd,fileuu C common /parablock/iam,nproc common /globsize/idimH(mzslp) common /procdat/ictxt,myrow,mycol common /contsize/icontold,nonzer common /matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest X ,NDIV common /allocblock/iflagallhbb,iflagallhbc,iflagallcf iflagallhbb=0 iflagallhbc=0 iflagallcf=0 c **** parallel **** C----------------------------------------------------------------------- IREAD = 5 IWRITE = 6 C ic = icolour/100 id = (icolour - 100*ic)/10 iu = (icolour - 100*ic - 10*id) ich0 = ichar('0') ic = ic + ich0 id = id + ich0 iu = iu + ich0 filecc = char(ic) filedd = char(id) fileuu = char(iu) C if(iam.eq.0)OPEN (UNIT=IWRITE, Z FILE='rout3r'//filecc//filedd//fileuu,STATUS='UNKNOWN', A FORM='FORMATTED') if(iam.ne.0)OPEN (UNIT=IWRITE, Z FILE='rout3rg'//filecc//filedd//fileuu,STATUS='UNKNOWN', A FORM='FORMATTED') C EX=.TRUE. c **** parallel **** inquire (file='dstg3',exist=ex) if (.not.ex) then write (iwrite,7777) istop=1 return endif 7777 format (/' dstg3 does not exist....stopping.') open(unit=iread,file='dstg3',status='old',form='formatted') c.......reading from sizeH.dat the size of the Hamiltonians call READ2R c write(0,*)'after read2r' if (istop.ne.0) return c **** parallel **** C----------------------------------------------------------------------- C C INITIALIZE COUNTER MSKIP: C MSKIP = 1 C C READ THE DATA SPECIFYING THE CASE C 10 CONTINUE c **** parallel **** CALL STG3RD C if (istop.ne.0) return c.... checking stop signal from process 0 NBUG10 = NBUG5 .NE. 0 .OR. NBUG6 .NE. 0 .OR. NBUG7 .NE. 0 .OR. A NBUG8 .NE. 0 IF (NBUG10 .AND. INAST.NE.0) WRITE (IWRITE,3000) C C c --- mcw: start --- nb0mcw=nb itestmcw=2*nb*max(nprow,npcol) if (itestmcw.gt.idimH(mskip)) then mcw=min(nprow,npcol) if(idimH(mskip).gt.mcw)then t=real(idimH(mskip)/mcw) else t=1 endif ipowmcw=int(log(t)/log(2.0d0))-1 if(ipowmcw.lt.0) ipowmcw=0 nb=2**ipowmcw mb=nb if(iam.eq.0) write(0,*) 'new NB: N,NB=',idimH(mskip), nb end if c --- mcw: end --- C C c...... reads overlap matrix now, if necessary. Here the size of the c (N+1)-Hamiltonian could be reduced. icontold=0 if (TOLER.ge.1.D-19) call HAMNEW0 c....... checking stop signal from process 0 if (istop.ne.0) return c.......initialize the BLACS context if (mskip.eq.1.and.nproctest.eq.0) call INITCTXT C write(0,*)iam,masteriam,mycol,myrow,ictxt C C C write(0,*)'myrow',myrow,'mycol',mycol,'iam',iam,'masteriam', C X masteriam c.......bails out if this process is not a part of this context if (myrow.eq.-1) return c.......calculate the size of the local Hamiltonian if(mskip.lt.ipwinit.or.idimH(mskip).eq.0)then irowsize=1 icolsize=1 allocate(hnp1(irowsize,icolsize)) goto 69 endif call HSIZE(irowsize,icolsize) c if(iam.eq.0)write(6,*)'hsize',icolour,mskip,more C C if(iam.eq.0) write(0,FMT='(3I6)')irowsize,icolsize,icolour c.......allocate the local Hamiltonian C allocate(HNP1(irowsize,icolsize),stat=ierr) if(ierr.ne.0)stop 'allocation fails - HNP1' C c.......initializate HNP1 do icol=1,icolsize do irow=1,irowsize HNP1(irow,icol)=0.0D0 enddo enddo c **** parallel **** C C C READ THE BASIC DATA AND THE HAMILTONIAN MATRIX FROM THE C FILE PRODUCED BY STG2 (NB: ON RETURN NOT2=MIN(NOT2,NRANG2)) C C 69 continue C C IF (IPLACE.NE.0) GOTO 70 C write(0,*)irowsize,icolsize,iam,masteriam,irsrc,icsrc c call comm_barrier() c call comm_finalize() c stop c **** parallel **** INQUIRE(file='STG2HCC000000',exist=esplit) IF(ESPLIT)THEN ioneone=0 !for safety CALL TAPERD2(HNP1,irowsize,icolsize) ELSE CALL TAPERD(HNP1,irowsize,icolsize) ENDIF c if(iam.eq.0)write(6,*)'taperd',icolour,mskip,more C if(iam.eq.0)write(0,*)'leaving taperd2',icolour call flush(6) c call comm_barrier() c call comm_finalize() c stop c.......checking stop signal from process 0 C C C if (istop.ne.0) return c **** parallel **** if((mskip.eq.ipwfinal).and.(icomplete_hdat.eq.1))then if(iam.eq.0)write(0,*)'full hdat write' more=0 endif c IF (IPLACE.NE.0.OR.MSKIP.LT.IPWINIT.or.idimH(mskip).eq.0)GOTO 40 CALL WRITOP IF (NBUG7.EQ.3) WRITE (IWRITE,3000) C C---- DIAGONALISE THE HAMILTONIAN (NB: RSCT RETURNS NOT1=MIN(=,NRANG2)) C c **** parallel **** C C CALL RSCT(HNP1,irowsize,icolsize) c write(0,*)'just called diag with more =',more,ifl c **** parallel **** IF (NBUG10) WRITE (IWRITE,3000) WRITE (IWRITE,3030) MSKIP C c if(npw_per_subworld.eq.1)call comm_barrier() !needed now? C C IF MORE.GT.0 THERE IS MORE DATA, INCREASE MSKIP BY 1, C IF MORE=0 THERE IS NO MORE DATA. C 40 CONTINUE deallocate(HNP1) !! parallel IF (MORE.EQ.0) GOTO 50 IF (NBUG10) WRITE (IWRITE,3010) 45 CONTINUE MSKIP = MSKIP + 1 C c --- mcw: start --- if(idimH(mskip).eq.0)go to 45 nb=nb0mcw mb=nb c --- mcw: end --- C C IF (MSKIP.LE.IPWFINAL) GOTO 10 C C---- PROCESS THE DIPOLE MATRIX ELEMENTS IF IPOLPH EQ 2 C 50 CONTINUE IF (IPOLPH.NE.2 .OR. IPLACE.GT.0 .OR. ITAPE1.EQ.0 X .OR. (IPLACE.LT.0.AND.INAST.NE.0) ) GO TO 60 INAST = MSKIP c CALL DMAT !! parallel !! DMAT not implemented 60 CONTINUE WRITE (IWRITE,3020) !! parallel 70 CONTINUE c CLOSE (IDISC1) !! HNP1 scratch not used by parallel RETURN C 3000 FORMAT (1X,71 ('*')) 3010 FORMAT (/' REPEAT STG3 WITH NEW DATA') 3020 FORMAT (/35X,'END OF STG3'/35X,11 ('-')) 3030 FORMAT (/' SYMMETRY NUMBER',I4,' COMPLETE.') END C********************************************************************** SUBROUTINE STG3RD C use param use basic2 use comm_interface IMPLICIT REAL*8 (A-H,O-Z) C C C----------------------------------------------------------------------- C C READS THE DATA SPECIFYING THE CASE C C----------------------------------------------------------------------- C C PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) PARAMETER (MXD1=MZNC2/MZNR2,MXD2=MZNR2/MZNC2,MXD3=MXD1+MXD2, A MXDM=MZNC2*MXD1/MXD3+MZNR2*MXD2/MXD3) C CHARACTER TITLE(18)*4,PARITY(0:1)*4,SPIN(0:9)*8,RAD*3 LOGICAL EX,EXA,EXH,ESPLIT,EXDIS C DIMENSION NTAPE(5) C COMMON /ADDE/EST(MZTAR),EPOLE,IPOLE,NEST,AST(MZTAR),MERGE(MZTAR) COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST,IONEONE COMMON /CUPINT/MNP1,NCONHP,NCHAN,N2HDAT,NTARSYM COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /INFORM/IREAD,IWRITE,IPUNCH COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /RCASES/NOT1,NOT2,NOT3,NPROG,NCFGP COMMON /RECOV/IPLACE COMMON /REDTOL/TOLER COMMON /SHIFT/ESHIFT,NSHIFT COMMON /NRBMNP/MNP2MX,IPRCENT,iprint COMMON /NRBSRT/ECORR(MZTAR),NCORR,ISORT COMMON /NRBHDX/EXH C NAMELIST/STG3A/NBUG5,NBUG6,NBUG7,NBUG8,NBUG9,IPOLPH,RAD,MNP2MX X ,ISORT,N2HDAT,NSINGLE,IPRCENT,iprint,ntarsym NAMELIST/STG3B/NBUT,NCONT,IDIAG,NAST,INAST,MAXC,TOLER,NSHIFT X ,ESHIFT,IPWINIT,IPWFINAL,NCORR,IONEONE C DATA NTAPE/1,2,3,4,5/,PARITY/'EVEN',' ODD'/,ZERO/0.0D0/, A SPIN/'(= 2J) ',' SINGLET',' DOUBLET',' TRIPLET',' QUARTET', B ' QUINTET',' SEXTET ',' SEPTET ',' OCTET ',' 2S+1>8 '/ c **** parallel **** include 'mpif.h' common /crashblock/istop c common /parablock/iam,nproc common /matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest X ,NDIV namelist/matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest A ,NDIV character*1 filec,filed,fileu character*1 filecc,filedd,fileuu c **** parallel **** c **** parallel **** C----------------------------------------------------------------------- C C IF ANY DIMENSION IS EXCEEDED WHEN READING FROM DATA, CALL RECOV2 C WITH IPLACE=0 TO TERMINATE THE PROGRAM C IPLACE = 0 C C IF MSKIP.GT.1, STG3 IS BEING REPEATED FOR NEW (N+1)-ELECTRON DATA C SO IT IS NOT NECESSARY TO READ IN THIS FIRST SET OF DATA AGAIN. C IF (MSKIP.GT.1) GOTO 10 C----------------------------------------------------------------------- MORE = 1 READ (IREAD,3120) TITLE C WRITE (IWRITE,3020) TITLE WRITE (IWRITE,3110) A MZCHF,MZLMX,MZLR2,MZMEG,MZKIL,MZNC2,MZNR2,MZSLP,MZTAR,MXMAT C C ICOPY = 0 ITOTAL = 0 IPUNCH=0 NBUG1=0 NBUG2=0 NBUG3=0 NBUG4=0 C IPOLPH=1 RAD='NO' NBUG5=0 NBUG6=0 NBUG7=0 NBUG8=0 NBUG9=0 C NSINGLE=0 MNP2MX=999999 !MAX MNP2, TO AID MEMORY SET ASIDE ISORT=0 N2HDAT=-1 !DEFAULT READ FROM FILE IPRCENT=100 !NON-PARTITIONED iprint=0 ntarsym=0 C READ(IREAD,STG3A) C IF(IPRCENT.GT.100)IPRCENT=100 IF(IPRCENT.LE.0)IPRCENT=60 !PARTITIONED IPOLPH = ABS(IPOLPH) IF(RAD.EQ.'YES')IPOLPH=2 c **** dario **** c options not (fully) implemented yet if (NBUG5.ne.0) then if(iam.eq.0)write(0,*) ' NBUG5.ne.0 not fully implemented' write(iwrite,*) ' NBUG5.ne.0 not fully implemented' endif if (NBUG6.ne.0) then if(iam.eq.0)write(0,*) ' NBUG6.ne.0 not implemented' write(iwrite,*) ' NBUG6.ne.0 not implemented' NBUG6=0 endif if (NBUG7.ne.0) then if(iam.eq.0)write(0,*) ' NBUG7.ne.0 not implemented' write(iwrite,*) ' NBUG7.ne.0 not implemented' NBUG7=0 endif if (NBUG8.ne.0) then if(iam.eq.0)write(0,*) ' NBUG8.ne.0 not implemented' write(iwrite,*) ' NBUG8.ne.0 not implemented' NBUG8=0 endif if (RAD.ne.'NO') then if(iam.eq.0)write(0,*) ' RAD.eq.YES not implemented' write(iwrite,*) ' RAD.eq.YES not implemented' RAD='NO' endif if (IPOLPH.ne.1) then if(iam.eq.0)write(0,*) ' IPOLPH.ne.1 not implemented' write(iwrite,*) ' IPOLPH.ne.1 not implemented' IPOLPH=1 endif c **** dario **** C C----------------------------------------------------------------------- C C Set the I/O numbers and open files. C C IREAD (5) .. input data .. dstg3 !! parallel C IWRITE (6) .. printed output .. rout3r C C IPUNCH .. NOT USED C C IDISC1 (11) .. scratch C IDISC2 (12) .. scratch (DA2) C IDISC3 (13) .. scratch C IDISC4 .. NOT USED C C ITAPE1 (1) .. RECUPD/STG2 dump (dipole) .. RECUPD/STG2D.DAT C if IPOLPH=2 C ITAPE2 (2) .. RECUPD/STG2 dump (hamiltonians) .. C RECUPH/STG2H.DAT/STG2HXXX.DAT C always used C ITAPE3 (3) .. STG3 dump (hamiltonians) .. H.DAT C always used C ITAPE4 (4) .. STG3 dump (dipole matrices) .. D00.DAT etc C if IPOLPH=2 C C----------------------------------------------------------------------- IWRITE = 6 IPUNCH = 0 C IDISC1 = 11 IDISC2 = 12 IDISC3 = 13 IDISC4 = 0 C ITAPE2 = 2 ITAPE3 = 3 IF (IPOLPH.EQ.2) THEN ITAPE1 = 1 ITAPE4 = 4 C FIX FOR DMAT1 TO AVOID DIMENSION ERROR AFTER LARGE CPU TIME USED. C SHOULD REALLY CODE MXDM=MAX(MZNC2,MZNR2,MZCHF) - NRB IF(MZCHF.GT.MXDM)CALL RECOV2('STG3RD','MZNC2 ',MZNC2,MZCHF+1) if (istop.ne.0) return ELSE ITAPE1 = 0 ITAPE4 = 0 ENDIF C c OPEN (UNIT=IDISC1,STATUS='SCRATCH',FORM='UNFORMATTED') !parallel C IF (ITAPE1.GT.0) THEN INQUIRE (FILE='RECUPD.DAT',EXIST=EX) IF (.NOT.EX) THEN INQUIRE (FILE='STG2D.DAT',EXIST=EX) IF (.NOT.EX) THEN WRITE (IWRITE,7771) 7771 FORMAT (/' Neither RECUPD.DAT nor STG2D.DAT exist....stopping.') istop=1 return ELSE OPEN (UNIT=ITAPE1,FILE='STG2D.DAT',STATUS='OLD', A FORM='UNFORMATTED') ENDIF ELSE OPEN (UNIT=ITAPE1,FILE='RECUPD.DAT',STATUS='OLD', A FORM='UNFORMATTED') ENDIF ENDIF INQUIRE (FILE='RECUPH000',EXIST=EX) IF (.NOT.EX) THEN INQUIRE (FILE='RECUPH.DAT',EXIST=EX) IF (.NOT.EX) THEN INQUIRE (FILE='STG2H.DAT',EXIST=EX) IF (.NOT.EX) THEN INQUIRE (FILE='STG2H000.DAT',EXIST=EX) IF(.NOT.EX) THEN WRITE (IWRITE,7772) 7772 FORMAT (/' Neither RECUPH.DAT nor STG2H.DAT/STG2H000.DAT exist.', A ' ...stopping.') istop = 3 return ELSE OPEN (UNIT=ITAPE2,FILE='STG2H000.DAT',STATUS='OLD', A FORM='UNFORMATTED') WRITE(IWRITE,3210) 3210 FORMAT(' OPENING FILE STG2H000.DAT'/) C INQUIRE (FILE='HDCORR000',EXIST=EXH) IF(EXH)OPEN(31,FILE='HDCORR000',STATUS='OLD', X FORM='FORMATTED') IF(EXH)WRITE(IWRITE,3212) 3212 FORMAT(' OPENING FILE HDCORR000'/) C INQUIRE (FILE='AMATU000.DAT',EXIST=EXA) IF(EXA)OPEN (UNIT=77,FILE='AMATU000.DAT', A STATUS='OLD',FORM='UNFORMATTED') IF(EXA)WRITE(IWRITE,3211) 3211 FORMAT(' OPENING FILE AMATU000.DAT'/) END IF ELSE OPEN (UNIT=ITAPE2,FILE='STG2H.DAT',STATUS='OLD', A FORM='UNFORMATTED') C INQUIRE (FILE='HDCORR.DAT',EXIST=EXH) IF(EXH)OPEN(31,FILE='HDCORR.DAT',STATUS='OLD', X FORM='FORMATTED') C INQUIRE (FILE='AMATU.DAT',EXIST=EXA) IF(EXA)OPEN (UNIT=77,FILE='AMATU.DAT',STATUS='OLD', A FORM='UNFORMATTED') ENDIF ELSE OPEN (UNIT=ITAPE2,FILE='RECUPH.DAT',STATUS='OLD', A FORM='UNFORMATTED') ENDIF ELSE OPEN (UNIT=ITAPE2,FILE='RECUPH000',STATUS='OLD', A FORM='UNFORMATTED') WRITE(IWRITE,3209) 3209 FORMAT(' OPENING FILE RECUPH000'/) END IF C C----------------------------------------------------------------------- WRITE (IWRITE,3000) WRITE (IWRITE,3010) A IREAD,IWRITE,IPUNCH, B IDISC1,IDISC2,IDISC3,IDISC4, C ITAPE1,ITAPE2,ITAPE3,ITAPE4 WRITE (IWRITE,3030) NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7, A NBUG8,NBUG9 WRITE (IWRITE,3050) ICOPY,ITOTAL,IPOLPH C IF (ITAPE1.GT.0) WRITE (IWRITE,3070) NTAPE(1),ITAPE1 WRITE (IWRITE,3070) NTAPE(2),ITAPE2 WRITE (IWRITE,3080) NTAPE(3),ITAPE3 IF (ITAPE4.GT.0) WRITE (IWRITE,3080) NTAPE(4),ITAPE4 C IF(IPRCENT.NE.100)WRITE(IWRITE,3085)IPRCENT C WRITE (IWRITE,3040) GOTO 20 C----------------------------------------------------------------------- 10 CONTINUE IF (INAST.EQ.0) GOTO 50 WRITE (IWRITE,3000) IF (INAST.GT.0) GOTO 40 C----------------------------------------------------------------------- C C READ AND WRITE BASIC DATA. C C NOT1 THE NUMBER OF CONTINUUM TERMS TO BE USED (DEFAULT=ALL) C IDIAG .EQ. 1 IF ITAPE2 CONTAINS STG2 HAMILTONIAN MATRICES, C AND THE DIAGONALIZATION DATA IS TO BE STORED ON ITAPE3 C NAST .EQ. 0 IF THE ENERGIES OF THE ATOMIC OR IONIC STATES ARE NOT C TO BE CHANGED *** C .NE. 0 IF THE ENERGIES ARE TO BE CHANGED. C THE VALUE OF NAST IS THEN THE TOTAL NUMBER OF C N-ELECTRON STATES (SEE BELOW FOR INPUT DETAILS). C .GT. 0 THEN IF EST(1).EQ.0, C EST(I)=OBSERVED ENERGIES RELATIVE TO THE GROUND, C IN RYDBERGS (EST(I).GT.0) OR /CM (EST(I).LT.0). C IF EST(1).NE.0, C EST(I)=OBSERVED ABSOLUTE ENERGIES IN A.U. C .LT. 0 MERGE GROUPS OF ENERGIES. C NCORR .GT. 0 READ NCORR TARGET INDEXES AND ENERGY CORRECTIONS C ECORR THAT WILL BE ADDED TO H_DIAG. IF NAST.GT.0 THEY FOLLOW C THOSE OBS. ENERGIES, AND ARE IN THE SAME UNITS, ELSE RYD. C ESHIFT .EQ. SHIFT OF ALL TARGET STATES IN RYDBERGS, CASE OF NAST .GT. 0 C EST(1)=0, TO MOVE THEM RELATIVE TO N+1 CFGS. C NSHIFT .GT.0 THEN READ FROM FILE 'eshift' THE SLPI, THE NUMBER OF C N+1-CFG EIGENVALUES TO BE SHIFTED, AND THEN THE INDEX C REALTIVE TO THE GROUND CONTINUUM OF THE N+1 EIGENVALUE C AND THE SHIFT (IN RY) OF THE EIGENVALUE - TWG C INAST .GE.0 IF STG3 IS TO BE RECYCLED WITH THE SAME BASIC DATA C FOR DIFFERENT (N+1)-ELECTRON SYMMETRIES (FOR ALL IF C .EQ.0, FOR INAST SPECIFIED SYMMETRIES IF .GT.0). C TOLER IS THE TOLERANCE PARAMETER FOR KEEPING/DISCARDING CERTAIN C LINEAR COMBINATIONS OF CONFIGURATIONS. C THE DEFAULT IS -1.0 (OFF) C IF YOU WANT TO KEEP ALL CFGS. SET THIS TO BE LARGER (E.G. 0.15) C IF YOU WANT NO (N+1)-CFGS., SET THIS TO BE 2. OR MORE - TWG & NRB C IPWINIT and C IPWFINAL USED TO SPECIFY AN INITIAL AND FINAL PARTIAL WAVE NUMBER C FROM THE LIST OF PARTIAL WAVES READ FROM STG2. THE CODE C WILL PROCESS ALL PARTIAL WAVES FROM IPWINIT TO IPWFINAL. C DEFAULT: ALL PARTIAL WAVES READ FROM STG2 ARE PROCESSED. C IF IPWINIT IS .GT. 1, THEN THE CODE WILL APPEND THE C NEWLY PROCESSED H.DAT FILE TO THE END OF THE EXISTING C H.DAT FILE. C----------------------------------------------------------------------- C 20 CONTINUE C NBUT=1 IDIAG=1 NAST=0 NCONT=999 MAXC=999 INAST=0 TOLER=-1.0D0 ESHIFT=0.0D0 NSHIFT=0 IPWINIT=1 IPWFINAL=MZSLP NCORR=0 IONEONE=0 C READ(IREAD,STG3B) c c IPWINIT= icolour*npw_per_subworld+1 IPWFINAL=(icolour+1)*npw_per_subworld C c READ(IREAD,STG3B) !surely we can't allow a user to redefine ipwinit,ipwfinal??? C if(iam.eq.0)write(0,21)'icolour= ',icolour,' ipwinit= ',ipwinit, X ' ipwfinal= ',ipwfinal,' npw_per_subworld= ',npw_per_subworld, X ' nproc=',nproc c 21 FORMAT(A9,I3,A10,I3,A11,I3,A19,I3,A7,I5) c c c **** dario **** c options not implemented yet if (IDIAG.ne.1) then if (iam.eq.0) write(0,*) ' IDIAG.ne.1 not implemented' write(iwrite,*) ' IDIAG.ne.1 not implemented' c IDIAG=1 istop=7 return endif if (INAST.ne.0) then if (iam.eq.0) write(0,*) ' INAST.ne.0 not implemented' write(iwrite,*) ' INAST.ne.0 not implemented' c INAST=0 istop=7 return endif c **** dario **** C ic = icolour/100 id = (icolour - 100*ic)/10 iu = (icolour - 100*ic - 10*id) ich0 = ichar('0') ic = ic + ich0 id = id + ich0 iu = iu + ich0 filecc = char(ic) filedd = char(id) fileuu = char(iu) C c IF (IPWINIT.EQ.1) THEN OPEN(UNIT=ITAPE3,FILE='H.DAT'//filecc//filedd//fileuu, + STATUS='UNKNOWN', + FORM='UNFORMATTED') c ELSE !why would we want to append now, when split, since 000 is skipped??? c OPEN(UNIT=ITAPE3,FILE='H.DAT'//filecc//filedd//fileuu, c + STATUS='UNKNOWN', cC + FORM='UNFORMATTED',ACCESS='APPEND') !f77 only c + FORM='UNFORMATTED',POSITION='APPEND') !f90 only c INAST=0 c ENDIF call comm_barrier() ESHIFT=ESHIFT/2.0D0 C IF(NSHIFT.GT.0) THEN OPEN(UNIT=81,FILE='ESHIFT.DAT',STATUS='OLD',FORM='FORMATTED') OPEN(UNIT=79,FILE='BBEIGEN.DAT',STATUS='OLD',FORM='FORMATTED') ENDIF C IF(TOLER.GE.1.D-19.AND..NOT.EXA)THEN WRITE(IWRITE,*) X 'FILE AMATU.DAT DOES NOT EXIST BUT TOLER.GT.0' if(iam.eq.0) x write(0,*)'FILE AMATU.DAT DOES NOT EXIST BUT TOLER.GT.0' istop=8 return ENDIF C c **** parallel **** itype=1 irsrc=0 icsrc=0 NB=32 MB=NB NPROW=-1 NPCOL=-1 nproctest = 0 !! test-run for checking dimensions NDIV=5 C read(iread,matrixdat) C C DISTRIBUTE.DAT override C INQUIRE(file='DISTRIBUTE.DAT',exist=exdis) IF(.NOT.exdis)then CONTINUE ELSE C OPEN(444,file='DISTRIBUTE.DAT',form='formatted',status='unknown') DO II=0,NPW-1 READ(444,*)IDUM,IDUM2,IDUM3,IDUM4 IF(ICOLOUR.eq.II)THEN NPROW=IDUM3 NPCOL=IDUM4 ENDIF ENDDO C CLOSE(444) ENDIF call MPI_Bcast(nprow,1,mpi_integer,0,new_comm,ier) call MPI_Bcast(npcol,1,mpi_integer,0,new_comm,ier) call mpi_barrier(new_comm,ierr) C C c........ dimension test mode: if (nproctest.ne.0) then nproc=nproctest t = nproc nprow = nint(sqrt(t)) npcol = nprow NOT1=MIN(NCONT,MAXC) NOT2=NOT1 return endif C IF(NPROW.LE.0.OR.NPCOL.LE.0)THEN T=NPROC NPROW=NINT(SQRT(T)) NPCOL=NPROW ENDIF if (nprow*npcol.ne.nproc) then if (iam.ne.0)write(0,*) + 'USAGE: nprow*npcol must be equal to # of processors' write(iwrite,*) + 'USAGE: nprow*npcol must be equal to # of processors' istop=1 return endif if (MB.ne.NB) MB=NB !FOR NOW c....... open timeXX.dat file (XX is the number of processors) if(iam.eq.0)then ic = nproc/100 id = (nproc - 100*ic)/10 iu = (nproc - 100*ic - 10*id) ich0 = ichar('0') ic = ic + ich0 id = id + ich0 iu = iu + ich0 filec = char(ic) filed = char(id) fileu = char(iu) if (ipwinit.eq.1) then open(unit=35,file='time'//filec//filed//fileu//'.dat', + status='unknown') write(35,2500) nproc,NB else open(unit=35,file='time'//filec//filed//fileu//'.dat', + status='unknown',position='append') write(35,2505) endif 2500 format(15x,'Timing'/3x,i3,' processors - (block size=:',i3,')'/) 2505 format(/5x,'New Run'/) endif c **** parallel **** NOT1=MIN(NCONT,MAXC) NOT2=NOT1 C WRITE (IWRITE,3100) NBUT,NOT1,NOT2,IDIAG,NAST,INAST C IF(EXH)write(IWRITE,*)'Using HDCORR to alter H diagonal' C C READ THE OBSERVED ENERGY OF THE ATOMIC OR IONIC TARGET STATES C NEST = NAST IF (NAST.EQ.0) GOTO 30 C C MERGE A GROUP OF ENERGIES TOGETHER C SEE SUBROUTINE TAPERD (GROUPING ENERGIES AROUND MEAN) C IF (NEST.LT.0) THEN READ(IREAD,*) (MERGE(I),I=1,-NEST) WRITE(IWRITE,*) -NEST, * ' target state groups, number of states in each group' WRITE(IWRITE,*) (MERGE(I),I=1,-NEST) GOTO 30 ENDIF C DO N=NAST+1,MZTAR EST(N)=ZERO ENDDO C READ (IREAD,*) (EST(N),N=1,NAST) C IF (EST(1).NE.ZERO) THEN WRITE (IWRITE,*) ' OBSERVED ENERGIES FOR EACH TARGET STATE (AU)' ELSE WRITE (IWRITE,*) ' OBSERVED ENERGIES RELATIVE TO GROUND STATE' IF (EST(NAST).GT.ZERO) WRITE (IWRITE,*) ' (RYD UNITS)' IF (EST(NAST).LT.ZERO) WRITE (IWRITE,*) ' (CM-1 UNITS)' ENDIF C WRITE (IWRITE,3090) (EST(N),N=1,NAST) C 30 CONTINUE C IF(NCORR.GT.0)THEN IF(NAST.EQ.0)THEN DO N=1,MZTAR EST(N)=ZERO ENDDO ENDIF DO N=1,MZTAR ECORR(N)=ZERO ENDDO DO N=1,NCORR READ (IREAD,*)I,ECORR(I) ENDDO ENDIF C----------------------------------------------------------------------- C C READ THE TOTAL ANGULAR MOMENTUM, SPIN, PARITY; C N.B. DO SET INAST=0 WHEN PROCESSING ALL SYMMETRIES! C IF (INAST.EQ.0) GOTO 50 C 40 CONTINUE C READ (IREAD,*) NSPN,LRGL,NPTY C LS = MIN(NSPN,9) WRITE (IWRITE,3060) LRGL,SPIN(LS),PARITY(NPTY) IF (NSPN.LT.0 .OR. NPTY.LT.0 .OR. NPTY.GT.1) then istop=1 !! parallel return endif IF (MSKIP.EQ.INAST) MORE = 0 C 50 CONTINUE C 3000 FORMAT (//35X,'SUBROUTINE STG3RD'/35X,17 ('-')) 3010 FORMAT ( A /' INPUT-OUTPUT CHANNEL NUMBERS' B /' ----------------------------' C//' IREAD (',I2,') .. input data .. dstg3' C /' IWRITE (',I2,') .. printed output .. rout3r' C /' IPUNCH (',I2,') .. NOT USED' C /' IDISC1 (',I2,') .. scratch' C /' IDISC2 (',I2,') .. scratch (DA2) ' C /' IDISC3 (',I2,') .. scratch' C /' IDISC4 (',I2,') .. NOT USED' C//' ITAPE1 (',I2,') .. RECUPD/STG2 dump (dipole) ', C '.. RECUPD/STG2D.DAT' C /' if IPOLPH=2' C//' ITAPE2 (',I2,') .. RECUPD/STG2 dump (hamiltonians) ', C '.. RECUPH/STG2H.DAT' C /' always used' C//' ITAPE3 (',I2,') .. STG3 dump (hamiltonians) ', C '.. H.DAT' C /' always used' C//' ITAPE4 (',I2,') .. STG3 dump (dipole matrices) ', C '.. D00 etc' C /' if IPOLPH=2' C ) 3020 FORMAT (//1X,72 ('-')//1X,18A4//1X,72 ('-')////15X, A ' SSSSSSSS TTTTTTTTTT GGGGGGGG 33333333 '/ B 15X, C 'SSSSSSSSSS TTTTTTTTTT GGGGGGGGGG 3333333333'/ D 15X, E 'SS TT GG GG 33 33'/ F 15X, G 'SS TT GG 33'/ H 15X, I 'SS TT GG 33'/ J 15X, K 'SSSSSSSSS TT GG 3333333 '/ L 15X, M ' SSSSSSSSS TT GG GGGG 3333333 '/ N 15X, O ' SS TT GG GGGG 33'/ P 15X, Q ' SS TT GG GG 33'/ R 15X, S ' SS TT GG GG 33 33'/ A 15X, B 'SSSSSSSSSS TT GGGGGGGGGG 33 33'/ C 15X, D ' SSSSSSSS TT GGGGGGGG 33333333 ') 3030 FORMAT (' NBUG PARAMETERS'/9I5) 3040 FORMAT ( A/' BASIC DATA' B/' ----------') 3050 FORMAT ( A/' ICOPY = ',I3 B/' ITOTAL = ',I3 C/' IPOLPH = ',I3) 3060 FORMAT (/' L =',I3,3X,A8,3X,A4/1X,24 ('-')) 3070 FORMAT (38X,' INPUT CHANNEL ITAPE',I1,' =',I5) 3080 FORMAT (38X,'OUTPUT CHANNEL ITAPE',I1,' =',I5) 3085 FORMAT (/'***** PARTITIONED R-MATRIX IN USE: ',I2, X'% OF E-SOLUTIONS IN USE.'/) 3090 FORMAT (5F16.6) 3100 FORMAT (/' NBUT =',I3,' NOT1 =',I3,' NOT2 =',I3,' IDIAG =',I3, A ' NAST =',I3,' INAST =',I3) 3110 FORMAT (//10X,'COMPILED FOR DIMENSIONS -'//15X, A 'MAX. NUMBER OF CHANNELS (NCHAN) MZCHF = ',I6/15X, C 'NO. OF MULTIPOLES IN POTENTIALS (LAMAX) MZLMX = ',I6/15X, E 'NO OF CONTINUUM ANGULAR MOM. (LRANG2) MZLR2 = ',I6/15X, G 'MEGAWORDS OF MEMORY AVAILABLE MZMEG = ',I6/15X, G 'KILOWORDS OF MEMORY AVAILABLE MZKIL = ',I6/15X, I 'N+1 ELECTRON CONFIGURATIONS MZNC2 = ',I6/15X, K 'CONTINUUM ORBITALS (NRANG2) MZNR2 = ',I6/15X, O 'NO OF (N+1) ELECTRON SYMMETRIES MZSLP = ',I6/15X, M 'TARGET STATES (NAST) MZTAR = ',I6/15X, P 'MAX. ORDER OF HAMILTONIAN MATRIX (MNP1) MXMAT = ',I6//) 3120 FORMAT (18A4) END C********************************************************************** SUBROUTINE TAPERD2(HNP1,irowsize,icolsize) use comm_interface C----------------------------------------------------------------------- C C C C----------------------------------------------------------------------- C C READS THE INPUT DATA FROM STG2 CORRESPONDING TO THE REQUIRED C VALUES OF LRGL, NSPN, NPTY, NCFGP. C C----------------------------------------------------------------------- C c use param use basic2 use boundH IMPLICIT REAL*8 (A-H,O-Z) INTEGER*8 JBIG,JJ,IJJ,IPOS,KP,KQ,KN,MO,MP,LBIG,NQ,MN2BIG,MNP1BIG INTEGER*8 :: L2PBIG(MZCHF) C PARAMETER (MXN21=MZNR2+1) PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) c PARAMETER (MXTRI=(MXMAT*(MXMAT+1))/2) c PARAMETER (MXD1=MZNC2/MZNR2,MXD2=MZNR2/MZNC2,MXD3=MXD1+MXD2, c A MXDM=MZNC2*MXD1/MXD3+MZNR2*MXD2/MXD3) c PARAMETER (MXDMAT=4*MXMAT*MZCHF+2*MXMAT*MXDM+2*MZNR2*MXDM, c A MXRSCT=MXTRI+MZCHF*MXMAT+MZNR2*MXMAT) c PARAMETER (MXM1=MXDMAT/MXRSCT,MXM2=MXRSCT/MXDMAT,MXM3=MXM1+MXM2, c A MXMATX=MXDMAT*MXM1/MXM3+MXRSCT*MXM2/MXM3+1) c PARAMETER (MXDUM2=MXMATX-MXTRI-MZNC2*MZNR2) C LOGICAL EX,EXA,EXBP,EXH,EXX CHARACTER PARITY(0:1)*4,SPIN(0:9)*8 CHARACTER*6 FL8(8),FO8(2) CHARACTER*1 IFNM(10),NUMA(0:9) CHARACTER*12 IFLNM CHARACTER*9 IFLBP CHARACTER*13 IFLNMCC CHARACTER*3 index1,index2,index3 C COMMON /ADDE/EST(MZTAR),EPOLE,IPOLE,NEST,AST(MZTAR),MERGE(MZTAR) COMMON /BASIC1/BSTO,RA,NELC,NZ,LRANG2,NRANG2,MAXNHF(MZLR2), A MAXNLG(MZLR2),MAXNC(MZLR2),LRANG1 COMMON /BASIN/EIGENS(MZNR2,MZLR2),ENDS(MXN21,MZLR2),DELTA,ETA COMMON /BUTT/COEFF(3,MZLR2),EK2MAX,EK2MIN,MAXNCB(MZLR2),NELCOR COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST,IONEONE COMMON /CUPINT/MNP1,NCONHP,NCHAN,N2HDAT,NTARSYM COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /INFORM/IREAD,IWRITE,IPUNCH c **** parallel **** c COMMON /MATRIX/HNP1(MXTRI),STORE(MZNC2,MZNR2),DUM2(MXDUM2) c **** parallel **** COMMON /MORDER/LAMAX,LAMBC,LAMCC,LAMIND,LAM COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /RCASES/NOT1,NOT2,NOT3,NPROG,NCFGP COMMON /RECOV/IPLACE COMMON /STATE/ENAT(MZTAR),LAT(MZTAR),ISAT(MZTAR),LPTY(MZTAR),NAST COMMON /SHIFT/ESHIFT,NSHIFT COMMON /NRBSRT/ECORR(MZTAR),NCORR,ISORT COMMON /NRBCOD/ICODE COMMON /NRBMNP/MNP2MX,IPRCENT,iprint COMMON /PART/TRACE COMMON /NRBHDX/EXH C DIMENSION JRELOP(3),NORDER(MZTAR) C SAVE IFL,nsym,KJI,KKJI C DATA IFL/1/,nsym/0/ DATA PARITY/'EVEN',' ODD'/,TWO/2.0D0/, A ZERO/0.0D0/,SPIN/'(= 2J) ',' SINGLET', A ' DOUBLET',' TRIPLET',' QUARTET',' QUINTET',' SEXTET ', B ' SEPTET ',' OCTET ',' 2S+1>8 '/ DATA FL8/' (8','(14X,7','(28X,6','(42X,5','(56X,4','(70X,3', A '(84X,2',' (98X,'/,FO8(2)/'F14.7)'/ DATA IFNM/'0','1','2','3','4','5','6','7','8','9'/ c **** parallel **** common /crashblock/istop include 'mpif.h' common /procdat/ictxt,myrow,mycol c common /parablock/iam,nproc dimension HNP1(irowsize,icolsize) allocatable store(:,:) allocatable haux(:) allocatable hcorr(:) COMMON /REDTOL/TOLER common /contsize/icontold,nonzer common /globsize/idimH(mzslp) character*1 filec,filed,fileu common /allocblock/iflagallhbb,iflagallhbc,iflagallcf c **** parallel **** C C----------------------------------------------------------------------- C C IF ANY DIMENSION IS EXCEEDED WHEN READING FROM DATA, CALL RECOV2 C WITH IPLACE=0 TO TERMINATE THE PROGRAM C NUMA(0)='0' NUMA(1)='1' NUMA(2)='2' NUMA(3)='3' NUMA(4)='4' NUMA(5)='5' NUMA(6)='6' NUMA(7)='7' NUMA(8)='8' NUMA(9)='9' C IPLACE = 0 WRITE (IWRITE,3000) ITAPE = ITAPE2 TRACE=ZERO !partitioned c nsym=nsym+1 c write(0,*)'mskip,nsym,ipolh=',mskip,nsym,IPOLPH inquire(file='sizeBP.dat',exist=exbp) if(.not.exbp.or.nsym.gt.1)then IF (MSKIP.GT.1 .AND. (INAST.EQ.0.OR. A ABS(IPOLPH).EQ.2)) GOTO 20 else IF (MSKIP.GT.1 .AND. (INAST.EQ.0.OR. A ABS(IPOLPH).EQ.2)) GOTO 8 endif c **** parallel **** C REWIND (77) !already called c **** parallel **** REWIND (ITAPE) C C READ THE BASIC DATA FROM FILE C 8 IF(N2HDAT.GT.0) THEN READ (ITAPE) NELC,NZ,LRANG1,LRANG2,NRANG2,LAMAX,ICODE,LAM,IZESP, A (JRELOP(I),I=1,3),MAXLD ELSE READ (ITAPE,ERR=9) NELC,NZ,LRANG1,LRANG2,NRANG2,LAMAX,ICODE,LAM, A IZESP,(JRELOP(I),I=1,3),MAXLD,N2HDT0 N2HDAT=N2HDT0 END IF 9 WRITE (IWRITE,3070) NELC,NZ,LRANG1, A LRANG2,NRANG2,LAMAX,LAM,JRELOP,IZESP IF (LRANG1.GT.MZLR2) CALL RECOV2('TAPERD',' MZLR2',MZLR2,LRANG1) IF (LRANG2.GT.MZLR2) CALL RECOV2('TAPERD',' MZLR2',MZLR2,LRANG2) IF (NRANG2.GT.MZNR2) CALL RECOV2('TAPERD',' MZNR2',MZNR2,NRANG2) IF (LAMAX.GT.MZLMX) CALL RECOV2('TAPERD',' MZLMX',MZLMX,LAMAX) if (istop.ne.0) return NOT2 = MIN(NRANG2,NOT2) !no longer used READ (ITAPE) (MAXNHF(L),L=1,LRANG1), (MAXNLG(L),L=1,LRANG1), A (MAXNC(L),L=1,LRANG1) WRITE (IWRITE,3050) (MAXNHF(L),L=1,LRANG1) WRITE (IWRITE,3060) (MAXNLG(L),L=1,LRANG1) DO 10 L = 1,LRANG2 READ (ITAPE) (EIGENS(N,L),N=1,NRANG2) READ (ITAPE) (ENDS(N,L),N=1,NRANG2) 10 CONTINUE READ (ITAPE) RA,BSTO,HINT,DELTA,ETA,NIX WRITE (IWRITE,3110) RA,BSTO,NIX C C IF(NIX.GT.0) THEN C READ(ITAPE) (IHX,I=1,NIX),(IRX,I=1,NIX) C IPTS=2*IRX C READ(ITAPE) (POVALU,I=1,IPTS) C DO 3 LP=1,LRANG1 C NBT=MAXNLG(LP)-LP+1 C IF(NBT.LE.0) GOTO 2 C DO 2 K=1,NBT C 2 READ(ITAPE) (POVALU,I=1,IPTS) C 3 CONTINUE C ENDIF if(exbp.AND.ifl.ne.1)then !Since RECUPDXXX complete, but don't want READ (ITAPE) !to overwrite ordered/observed energies READ (ITAPE) READ (ITAPE) goto 20 endif READ (ITAPE) ((COEFF(I,L),I=1,3),L=1,LRANG2) READ (ITAPE) NAST WRITE(IWRITE,*)' NAST = ',NAST NAST = ABS(NAST) IF (NAST.GT.MZTAR) CALL RECOV2('TAPERD',' MZTAR',MZTAR,NAST) if (istop.ne.0) return READ (ITAPE) (ENAT(I),I=1,NAST), (LAT(I),I=1,NAST), A (ISAT(I),I=1,NAST), (LPTY(I),I=1,NAST) C C THE FOLLOWING DATA DEPENDS ON LRGL, NSPN ETC. C C Initialise partial wave index (index2) C KKJI=-1 JI=-1 C C 20 CONTINUE C C KKJI=KKJI+1 KKJI1=KKJI/100 KKJI2=(KKJI-100*KKJI1)/10 KKJI3=KKJI-100*KKJI1-10*KKJI2 C index2=NUMA(KKJI1)//NUMA(KKJI2)//NUMA(KKJI3) C DO JI=0,NTARSYM-1 JI1=JI/100 JI2=(JI-100*JI1)/10 JI3=JI-100*JI1-10*JI2 C index1=NUMA(JI1)//NUMA(JI2)//NUMA(JI3) c INQUIRE(file='STG2HCC'//index1//index2,exist=exx) C if(.not.exx)then call mpi_barrier(mpi_comm_world,ierr) call mpi_finalize(ierr) stop endif OPEN(88,file='STG2HCC'//index1//index2,form='unformatted', X status='unknown') c if(iam.eq.0)write(0,*)'opening STG2HCC'//index1//index2 C C C C ******************** test for empty Hamiltonians ******* C C symmetry could be empty : test C READ(88,IOSTAT=IERROR)LRGL222 C if(ierror.eq.0)goto 22 C if(iam.eq.0)write(0,*)'closing STG2HCC'//index1//index2 CLOSE(88) C ENDDO C if(iam.eq.0)write(0,*)'no coupled channels for this symm' goto 20 C C 22 CONTINUE ierror=0 BACKSPACE(88) C C KJI=0 !(index over HCC blocks ) C C C ********************************************************* C READ (88,END=230) LRGL2,NSPN2,NPTY2,NCFGP2,IPOL2 !note WRITE (IWRITE,3130) LRGL2,NSPN2,NPTY2,NCFGP2 C C C IF((IFL.GT.1).AND.(IONEONE.eq.1))THEN WRITE(IWRITE,*)' ' WRITE(IWRITE,*)'SKIPPING EARLY STG2HXXX.DATS IF EARLY SYMMETRIES' WRITE(IWRITE,*)'ARE NOT REQUIRED and WE HAVE 1 SYM PER FILE' WRITE(IWRITE,*)' ' MSKIP=MSKIP+1 GOTO 210 ENDIF C C C C READ (88) MNP1,NCONHP,NCHAN !note c if(iam.eq.0)write(*,*)mnp1,nconhp,nchan IF (NCHAN.GT.MZCHF) CALL RECOV2('TAPERD',' MZCHF',MZCHF,NCHAN) IF (MNP1.GT.MXMAT) CALL RECOV2('TAPERD','MXMAT ',MXMAT,MNP1) IF (NCFGP2.GT.MZNC2) CALL RECOV2('TAPERD',' MZNC2',MZNC2,NCFGP2) if (istop.ne.0) return READ (88) (NCONAT(I),I=1,NAST) !note IF(NSPN2.NE.0)READ (88) (L2P(I),I=1,NCHAN) !note IF(NSPN2.EQ.0)READ (88) (L2P(I),I=1,NCHAN),(KJ(I),I=1,NCHAN) ! NRB READ (88) MORE2 !note c if(iam.eq.0)write(*,*)more2 c call flush_(iwrite) C C c if HDCORR.DAT exists, subtract corrections from diagonal c IF(EXH.AND.NCFGP2.GT.0)THEN READ(31,*)LRGL,NSPN,NPTY,mdiag IF(LRGL.NE.LRGL2.OR.NSPN.NE.NSPN2.OR.NPTY.NE.NPTY2)THEN WRITE(IWRITE,*) 'MIS-MATCH ON HDCORR.DAT FOR L S P:', X LRGL2,NSPN2,NPTY2,LRGL,NSPN,NPTY STOP 'MIS-MATCH ON HDCORR.DAT FOR L S P' ENDIF IF(MNP1.NE.mdiag)THEN WRITE(IWRITE,*) 'MIS-MATCH ON HDCORR.DAT FOR MNP1:', X MNP1,mdiag STOP 'MIS-MATCH ON HDCORR.DAT FOR MNP1' ENDIF c **** parallel **** allocate(hcorr(mnp1),stat=ierr) if(ierr.ne.0)stop'program could not allocate memory for hcorr' c **** parallel **** do j=1,MNP1 read(31,*)jjj,hcorr(j) enddo ENDIF C c **** parallel **** allocate(store(nrang2,nrang2),stat=ierr) if(ierr.ne.0)stop'program could not allocate memory for store' c **** parallel **** C C READ THE UPPER TRIANGLE OF THE HAMILTONIAN MATRIX INTO HNP1 C do M=1,NCHAN L2PBIG(M)=L2P(M) enddo C MDIAGCOUNT=1 !partitioned ICOUNTA=0 !partitioned mdiag=0 !hcorr C c if(iam.eq.0)write(*,*)'starting CC read',icode MN2 = 2*MNP1 MN2BIG=MN2 MNP1BIG=MNP1 C C LT=0 C if(npw_per_subworld.eq.1.and.mskip.ne.(icolour+1))then goto 111 endif if(mskip.lt.ipwinit)goto 111 C DO 110 M = 1,NCHAN MP = (M* (M+1))/2 MO = MP - M + 1 NQ = (M-1)*NRANG2 NBT = 1 DO 100 LBIG = MO,MP C C ---> READ CONTINUUM-CONTINUUM BLOCKS INTO STORE IN TURN IF (ICODE.NE.21) GOTO 50 IF (LBIG.EQ.MP) GOTO 40 IF (L2PBIG(M).EQ.L2PBIG(LBIG-MO+1)) GOTO 50 40 CONTINUE goto 48 C C C 45 CONTINUE if(KJI.eq.(ntarsym-1))goto 111 CLOSE(88) C C C KJI=KJI+1 KJI1=KJI/100 KJI2=(KJI-100*KJI1)/10 KJI3=KJI-100*KJI1-10*KJI2 index1=NUMA(KJI1)//NUMA(KJI2)//NUMA(KJI3) C C C C C C write(0,*)KJI,KJI1,KJI2,KJI3 OPEN(88,file='STG2HCC'//index1//index2, X form='unformatted',status='unknown') C C C if(iam.eq.0)then c write(0,*)'opening STG2HCC'//index1//index2 endif C C C 48 CONTINUE C READ(88,IOSTAT=IERROR)LRGL222 C IF(IERROR.EQ.0)THEN BACKSPACE(88) ELSE GOTO 45 ENDIF C if(iam.eq.0)then READ (88,END=45) ((STORE(I,J),I=1,NRANG2),J=1,NRANG2) else READ (88,END=45) endif C if(iam.eq.0)write(1000,*)((STORE(I,J),I=1,NRANG2),J=1,NRANG2) IF (LBIG.GE.MP) GOTO 70 GOTO 60 C C 50 CONTINUE C READ(88,IOSTAT=IERROR)LRGL222 C IF(IERROR.EQ.0)THEN BACKSPACE(88) ELSE GOTO 45 ENDIF C if(iam.eq.0)then READ (88,END=45) ((STORE(J,I),I=1,NRANG2),J=1,NRANG2) else READ (88,END=45) endif C if(iam.eq.0)write(1000,*)((STORE(I,J),I=1,NRANG2),J=1,NRANG2) 60 CONTINUE KP = (LBIG-MO)*NRANG2 KN = ((MN2BIG-KP+1)*KP)/2 + NQ - KP C C ---> TRANSFER CONTINUUM-CONTINUUM BLOCKS INTO HNP1 FROM STORE 70 CONTINUE N222=NRANG2*NRANG2 call MPI_BCAST(STORE,N222,MPI_REAL8,0,NEW_COMM,ierr) C C SKIP IF THIS SYMMETRY IS NOT NEEDED IF (MSKIP.LT.IPWINIT) GO TO 100 C DO 90 J = 1,NRANG2 IF (LBIG.NE.MP) THEN JBIG=J JJ = (JBIG-1)* (MNP1BIG-KP) - ((JBIG-1)*JBIG)/2 + KN C ELSE NBT = J JBIG=J JJ = ((NQ+JBIG-1)* (MN2BIG-NQ-JBIG+2))/2 - JBIG + 1 ENDIF C DO 80 I = NBT,NRANG2 c **** parallel **** c...........in serial stg3: HNP1(I+JJ) = STORE(I,J) c...........here I+JJ global position is mapping into local (irow,icol) if (toler.lt.1.d-19.or.nonzer.eq.0) then ipos=I+JJ else c...............transform the index I+JJ for size MNP1 to c...............the index ipos for size idimH(mskip) IJJ=I+JJ c call ntoij(mnp1,IJJ,ig,jg) c call ijton(idimH(mskip),ipos,ig,jg) call ntoijandback(mnp1,IJJ,ipos) endif C C Locate diagonal both for HCC contribution to TRACE & hcorr C if(IPRCENT.ne.100.or.EXH)then if(ipos.eq.mdiagcount)then if(EXH.AND.NCFGP2.GT.0)then mdiag=mdiag+1 store(i,j)=store(i,j)-hcorr(mdiag) endif if(IPRCENT.ne.100)TRACE=TRACE+store(i,j) !Partitioned mdiagcount=mdiagcount+(mnp1-icounta) icounta=icounta+1 endif endif c call locH(ipos,irow,icol,ipr,ipc) if (ipr.eq.myrow.and.ipc.eq.mycol) + HNP1(irow,icol) = store(i,j) C c **** parallel **** 80 CONTINUE 90 CONTINUE 100 CONTINUE 110 CONTINUE c 111 CONTINUE if(allocated(store))deallocate(store) C close(88) C call mpi_barrier(mpi_comm_world,ierr) C C C if(ncfgp2.ne.0)then OPEN(88,file='STG2HBC'//NUMA(0)//NUMA(0)//NUMA(0)//index2, X form='unformatted',status='unknown') endif C C C if(iam.eq.0)then c write(0,*)'opening STG2HBC'//NUMA(0)//NUMA(0)//NUMA(0)//index2 endif C c if(iam.eq.0)write(*,*)'ending CC read' C C READ BOUND-CONTINUUM BLOCKS INTO STORE C IF (NCFGP2.EQ.0) GOTO 160 IF (NCFGP2.GT.MZNC2) CALL RECOV2('TAPERD',' MZNC2',MZNC2,NCFGP2) if (istop.ne.0) return c **** parallel **** if (toler.ge.1.d-19.and.nonzer.gt.0.or.nshift.gt.0) then if (iflagallhbc.ne.0) then deallocate(hbc) iflagallhbc = 0 endif if (iflagallhbb.ne.0) then deallocate(hbb) iflagallhbb = 0 endif if(toler.ge.1.d-19.and.nonzer.gt.0)then ncont=mnp1-ncfgp2 allocate(hbc(ncfgp2,ncont),stat=ierr) if(ierr.ne.0)stop'allocation fails for hbc' iflagallhbc = 1 endif allocate(hbb(ncfgp2,ncfgp2),stat=ierr) if(ierr.ne.0)stop'allocation fails for hbb' iflagallhbb = 1 endif allocate(store(ncfgp2,nrang2),stat=ierr) if(ierr.ne.0) stop 'allocation fails for store (b)' c **** parallel **** c if(iam.eq.0)write(*,*)'starting BC read' C if(npw_per_subworld.eq.1.and.mskip.ne.(icolour+1))then goto 141 endif if(mskip.lt.ipwinit)goto 141 c DO 140 M = 1,NCHAN READ (88) ((STORE(I,J),I=1,NCFGP2),J=1,NRANG2) N222=NCFGP2*NRANG2 call MPI_BCAST(STORE,N222,MPI_REAL8,0,NEW_COMM,ierr) C if(iam.eq.0)write(2000,*)((STORE(I,J),I=1,NCFGP2),J=1,NRANG2) C C SKIP IF THIS SYMMETRY IS NOT NEEDED IF (MSKIP.LT.IPWINIT) GO TO 140 C KP = (M-1)*NRANG2 KQ = (KP* (MN2BIG-KP+1))/2 + NCONHP DO 130 J = 1,NRANG2 JBIG=J JJ = (JBIG-1)*MNP1BIG - (JBIG* (JBIG-1))/2 - JBIG*KP + KQ DO 120 I = 1,NCFGP2 c **** parallel **** c........in serial stg3: HNP1(I+JJ) = STORE(I,J) if (toler.lt.1.d-19.or.nonzer.eq.0) then c............here I+JJ global position is mapping into local (irow,icol) IJJ=I+JJ call locH(IJJ,irow,icol,ipr,ipc) if (ipr.eq.myrow.and.ipc.eq.mycol) + HNP1(irow,icol) = store(i,j) else c...............put in HBC for further transformation in HAMNEW IJJ=I+JJ call ntoij(mnp1,IJJ,irow,icol) icol=icol-icontold c................HBC is the transverse of HCB HBC(icol,irow) = store(i,j) endif c **** parallel **** 120 CONTINUE 130 CONTINUE 140 CONTINUE c 141 CONTINUE c deallocate(store) allocate(haux(ncfgp2),stat=ierr) if(ierr.ne.0)stop'allocation fails for haux' C c if(iam.eq.0)write(*,*)'ending BC read' C C READ BOUND-BOUND ELEMENTS INTO HNP1 C CLOSE(88) C C call mpi_barrier(mpi_comm_world,ierr) C OPEN(88,file='STG2HBB'//NUMA(0)//NUMA(0)//NUMA(0)//index2, X form='unformatted',status='unknown') if(iam.eq.0)then c write(0,*)'opening STG2HBB'//NUMA(0)//NUMA(0)//NUMA(0)//index2 endif c if(iam.eq.0)write(*,*)'starting BB read' c if(mdiag.ne.nconhp)stop 'mdiag.ne.nconhp' C MP = (NCONHP* (MN2BIG-NCONHP+1))/2 C if(npw_per_subworld.eq.1.and.mskip.ne.(icolour+1))then goto 151 endif if(mskip.lt.ipwinit)goto 151 C DO 150 I = 1,NCFGP2 MO = MP + 1 MP = MO + NCFGP2 - I c **** parallel **** c........in serial stg3: READ (ITAPE) (HNP1(J),J=MO,MP) c........here an auxiliar haux array is created which has hnp1(J) J=mo,mp c........then, the J global position is mapping into local (irow,icol) c read (88) (haux(j),j=1,mp-mo+1) N222=mp-mo+1 call MPI_BCAST(haux,N222,MPI_REAL8,0,NEW_COMM,ierr) C if(iam.eq.0)write(3000,*)(haux(j),j=1,mp-mo+1) c c........skip if this symmetry is not needed if (mskip.lt.ipwinit) go to 150 c jj=1 C Locate diagonal both for HBB contribution to TRACE & hcorr if(EXH.AND.NCFGP2.GT.0)haux(jj)=haux(jj)-hcorr(nconhp+i) if(IPRCENT.ne.100)TRACE=TRACE+haux(jj) !Partitioned c C do 149 j=mo,mp do 149 jbig=mo,mp c if ((toler.lt.1.d-19.or.nonzer.eq.0).and.nshift.le.0) then c call locH(j,irow,icol,ipr,ipc) call locH(JBIG,irow,icol,ipr,ipc) if (ipr.eq.myrow.and.ipc.eq.mycol) + HNP1(irow,icol) = haux(jj) else c.............put in HBB for further transformation in HAMNEW(2) c call ntoij(mnp1,J,irow,icol) call ntoij(mnp1,jbig,irow,icol) irow=irow-icontold icol=icol-icontold HBB(irow,icol) = haux(jj) HBB(icol,irow) = haux(jj) !symmetrize endif jj=jj+1 149 continue c **** parallel **** 150 CONTINUE c 151 CONTINUE c deallocate(haux) if(EXH.AND.NCFGP2.GT.0)deallocate(hcorr) C c if(iam.eq.0)write(*,*)'ending BB read' C C READ THE ASYMPTOTIC COEFFICIENTS (ARRAY CF EXPANDED IN WRITOP) C 160 CONTINUE C C The following covers the case of no N+1's but just a CF read C if(ncfgp2.eq.0)then OPEN(88,file='STG2HBB'//NUMA(0)//NUMA(0)//NUMA(0)//index2, X form='unformatted',status='unknown') endif C C C call mpi_barrier(mpi_comm_world,ierr) IF (LAMAX.EQ.0) GOTO 180 c c **** parallel **** if (iflagallcf.ne.0) then deallocate(cf) iflagallcf = 0 endif C C C if(mskip.eq.(icolour+1))then allocate(cf(nchan,nchan,lamax),stat=ierr) else allocate(cf(1,1,1)) endif C C C if(ierr.ne.0)stop'allocation fails for cf' iflagallcf = 1 c **** parallel **** c c if(mskip.eq.(icolour+1))then DO 170 I = 1,NCHAN READ (88) ((CF(I,J,K),K=1,LAMAX),J=I,NCHAN) 170 CONTINUE endif C C CHECK THAT LRGL,NSPN,NPTY,NCFGP HAVE THE REQUIRED VALUES. C IF NOT, READ IN MORE DATA FROM FILE OR USER. C C*******************************SPLIT*********************** CPB C C goto 230 C C 180 CONTINUE IF (MORE2.EQ.0.and..not.exbp) THEN MORE = 0 END IF IF (INAST.NE.0) GOTO 190 LRGL = LRGL2 NSPN = NSPN2 NPTY = NPTY2 190 CONTINUE c IF (LRGL.EQ.LRGL2 .AND. NSPN.EQ.NSPN2 .AND. A NPTY.EQ.NPTY2) THEN if(exbp.and.more2.eq.0)goto 210 GOTO 230 ENDIF c IF (IPOLPH.LE.1) GOTO 200 WRITE (IWRITE,3140) MORE = 0 GOTO 220 C 200 CONTINUE IF (MORE2.GE.1) THEN if (iflagallcf.ne.0) then deallocate(cf) iflagallcf = 0 endif if(exbp)go to 8 GOTO 20 ENDIF 210 CONTINUE C IF(N2HDAT.EQ.IFL) THEN !NO MORE FILES IF(INAST.EQ.0) THEN MORE=0 ELSE WRITE (IWRITE,3020) IF(MSKIP.LT.INAST) MORE=1 END IF ELSE C C OPEN NEW FILE STG2HXXX.DAT ich0 = ichar('0') ic = ifl/100 id = (ifl - 100*ic)/10 iu = (ifl - 100*ic - 10*id) ic = ic + ich0 id = id + ich0 iu = iu + ich0 filec = char(ic) filed = char(id) fileu = char(iu) c nsym=0 c if(.not.exbp)then c IFLNMCC='STG2HCC000'//filec//filed//fileu INQUIRE (FILE=IFLNM,EXIST=EX) IF(.NOT.EX) THEN WRITE(IWRITE,3200) IFLNM istop=4 return ELSE CLOSE(ITAPE) OPEN (UNIT=ITAPE,FILE=IFLNM,STATUS='OLD',FORM='UNFORMATTED') WRITE(IWRITE,3210) IFLNM END IF IFLNM='HDCORR'//filec//filed//fileu INQUIRE (FILE=IFLNM,EXIST=EXH) IF(EXH)THEN CLOSE(31) OPEN (31,FILE=IFLNM,STATUS='OLD',FORM='FORMATTED') WRITE(IWRITE,3210) IFLNM ENDIF IFL=IFL+1 GO TO 20 c else c IFLBP='RECUPH'//filec//filed//fileu INQUIRE (FILE=IFLBP,EXIST=EX) IF(.NOT.EX) THEN WRITE(IWRITE,3205) WRITE(IWRITE,3200) IFLBP WRITE(IWRITE,3205) WRITE(IWRITE,3220) WRITE(IWRITE,3205) MORE=0 ELSE CLOSE(ITAPE) OPEN (UNIT=ITAPE,FILE=IFLBP,STATUS='OLD',FORM='UNFORMATTED') WRITE(IWRITE,3205) WRITE(IWRITE,3210) IFLBP WRITE(IWRITE,3205) MORE=1 END IF IFL=IFL+1 c IPLACE = -1 GOTO 230 c endif c ENDIF c 220 CONTINUE IPLACE = -1 GOTO 390 230 CONTINUE C c if(iam.eq.0)write(0,*)'soundoff',icolour c call flush(0) C C C if(iam.eq.0)write(0,*)ifl,nsym,mskip,more nsym=0 IFL=IFL+1 NCFGP=NCFGP2 LRGL=LRGL2 NPTY=NPTY2 NSPN=NSPN2 CLOSE(88) C C JI=KKJI+1 JI1=JI/100 JI2=(JI-100*JI1)/10 JI3=JI-100*JI1-10*JI2 C index3=NUMA(JI1)//NUMA(JI2)//NUMA(JI3) C c INQUIRE(file='STG2HCC000'//index3,exist=ex) if(.not.ex)then MORE=0 if(iam.eq.0)write(0,*)'LAST SYMMETRY' else MORE=1 endif C c goto 235 C C SKIP IF THIS SYMMETRY IS NOT NEEDED c IF (MSKIP.LT.IPWINIT) GO TO 235 C C C TRANSFORM THE HAMILTONIAN TO REDUCE NUMBER OF N+1-CFGS - TWG 13/4/95. C c IF(NCFGP2.GT.0.and.TOLER.ge.1.D-19.and.nonzer.gt.0)THEN c call HAMNEW(HNP1,irowsize,icolsize, c + MNP1,NCFGP2,LRGL2,NSPN2,NPTY2) c MN2=2*MNP1 c ENDIF c NCFGP = NCFGP2 C C SHIFT THE N+1 CFGS - TWG 18/11/95. C c IF(NCFGP2.GT.0.and.NSHIFT.gt.0) THEN !parallel c CALL HAMNEW2(HNP1,irowsize,icolsize, c + MNP1,NCFGP2,LRGL2,NSPN2,NPTY2) c ENDIF c NCFGP = NCFGP2 !not redefined c 235 CONTINUE c **** parallel **** if(iflagallhbc.ne.0)then deallocate(hbc) iflagallhbc = 0 endif if(iflagallhbb.ne.0)then deallocate(hbb) iflagallhbb = 0 endif c **** parallel **** C C WRITE OUT THE TARGET STATES. C IF (MSKIP.GT.1) GOTO 250 WRITE (IWRITE,3030) CALL ORDER(ENAT,NORDER,NAST,1) IF (NCORR+ABS(NEST).EQ.0) THEN GROUND = ENAT(NORDER(1)) DO 240 II = 1,NAST I = NORDER(II) T = (ENAT(I)-GROUND)*2 IS = MIN(ISAT(I),9) WRITE (IWRITE,3040) SPIN(IS),LAT(I), A PARITY(LPTY(I)),NCONAT(I),ENAT(I),ZERO,T,II 240 CONTINUE ENDIF C C IF NEST=NAST, ADD ON THE EXTRA ENERGY TO EACH STATE AND TO THE C ASSOCIATED DIAGONAL ELEMENTS OF THE HAMILTONIAN MATRIX. C C EXTENSION '90JUL1/5: OBS.ENERGY INPUT IN RYD ABOVE 1 OR IN CM**-1 C EST(I)=0 ENSURES UNCORRECTED TERM ENERGY IF NO OBS.DATA AVAILABLE C 250 CONTINUE IF (NCORR+ABS(NEST).EQ.0) GOTO 310 E0 = ENAT(1) IF (MSKIP.LE.1) THEN C C IF NEST .LT. 0 GROUP ENERGIES AROUND ARITHMETIC MEANS C IF (NEST.LT.0) THEN DO II = 1,NAST EST(II) = ENAT(II) ENDDO ISTART = 0 DO 258 I = 1,-NEST IF(ISTART.GE.NAST.OR.MERGE(I).LE.0) GOTO 258 IF(MERGE(I).EQ.1) GOTO 256 EMEAN = 0.0D0 IEND = MIN(ISTART+MERGE(I),NAST) DO II=ISTART+1,IEND EMEAN = EMEAN+ENAT(II) ENDDO EMEAN = EMEAN / (IEND-ISTART) DO II=ISTART+1,IEND EST(II) = EMEAN ENDDO 256 ISTART = ISTART + MERGE(I) 258 CONTINUE C IF (ISTART.NE.NAST) THEN WRITE(IWRITE,*) '**** WARNING **** ', A 'Target State Energy MERGE only treats', B ISTART,' states' ENDIF ENDIF C C IF NCORR,ABS(NEST) .GT. 0 CORRECT ENERGIES C IF(EST(1).EQ.ZERO)E0=E0+ESHIFT DO 270 II = 1,NAST I = NORDER(II) I0=I IF(ISORT.GT.0)I0=II C AST(I) = ZERO IF(NEST.NE.0)AST(I) = EST(I0) - ENAT(I) IF(NCORR.GT.0)AST(I)=AST(I)+ECORR(I0) C IF (EST(1).EQ.ZERO) THEN AST(I) = ZERO IF(NEST.NE.0)AST(I) = EST(I0) IF(NCORR.GT.0)THEN IF(ECORR(I0).NE.0)THEN AST(I)=AST(I)+ECORR(I0) IF(EST(I0).EQ.0)AST(I)=AST(I)+(ENAT(I)-E0)*TWO ENDIF ENDIF AST(I)=AST(I)/TWO C IF (AST(I).EQ.ZERO)THEN ENAT(I) = ENAT(I) + ESHIFT GOTO 260 ENDIF C IF (EST(I0).LT.ZERO) AST(I) = -AST(I)/109737.3D0 AST(I) = AST(I) - ENAT(I) + E0 ENDIF C ENAT(I) = ENAT(I) + AST(I) C 260 CONTINUE T = (ENAT(I)-ENAT(1))*TWO IS = MIN(ISAT(I),9) WRITE(IWRITE,3040)SPIN(IS),LAT(I),PARITY(LPTY(I)),NCONAT(I), A ENAT(I),AST(I),T,II 270 CONTINUE ELSE IF(INAST.GT.0)THEN C ENAT HAS BEEN RE-READ AND SO MUST BE RE-ADJUSTED : NRB 4/11/94 CALL ORDER(ENAT,NORDER,NAST,1) DO 275 II=1,NAST I=NORDER(II) ENAT(I)=ENAT(I)+AST(I) 275 CONTINUE ENDIF C END MODS ENDIF C C C LOOP OVER THE CHANNELS ASSOCIATED WITH EACH STATE IN THE H-MATRIX C C SKIP IF THIS SYMMETRY IS NOT NEEDED IF (MSKIP.LT.IPWINIT) GO TO 390 IF (NEST.EQ.0) GOTO 310 C MN2BIG=MN2 !CPB oct 2005 (as MN2 has been redefined) C MP = 0 DO 300 I = 1,NAST IF (NCONAT(I).LE.0) GOTO 300 MO = MP + 1 MP = MP + NCONAT(I) DO 290 M = MO,MP NQ = (M-1)*NRANG2 DO 280 J = 1,NRANG2 JBIG=J JJ = ((NQ+JBIG-1)* (MN2BIG-NQ-JBIG+2))/2 + 1 c **** parallel **** c............in serial stg3: HNP1(JJ) = HNP1(JJ) + AST(I) call locH(JJ,irow,icol,ipr,ipc) if (ipr.eq.myrow.and.ipc.eq.mycol) + HNP1(irow,icol) = HNP1(irow,icol) + AST(I) c **** parallel **** 280 CONTINUE 290 CONTINUE 300 CONTINUE C C WRITE THE FULL HAMILTONIAN MATRIX TO SCRATCH DISC C 310 CONTINUE c **** parallel **** c all this chunk is commented - need to be implemented in the future c IF (NOT1.GE.NRANG2) GOTO 330 c REWIND IDISC1 c DO 320 M = 1,NCHAN c MO = ((M-1)*NRANG2* (MN2- (M-1)*NRANG2+1))/2 + 1 c MP = (M*NRANG2* (MN2-M*NRANG2+1))/2 c WRITE (IDISC1) (HNP1(I),I=MO,MP) c 320 CONTINUE c MO = (NCONHP* (MN2-NCONHP+1))/2 + 1 c MP = (MNP1* (MNP1+1))/2 c WRITE (IDISC1) (HNP1(I),I=MO,MP) C C WRITE OUT THE HAMILTONIAN MATRIX IF NBUG6=1 C c 330 CONTINUE WRITE (IWRITE,3090) MNP1,NCHAN c IF (NBUG6.NE.1) GOTO 370 c M = 0 c 340 CONTINUE c K = M c L = M + 1 c M = M + 8 c JJ = MIN(M,MNP1) c FO8(1) = FL8(1) c DO 360 I = 1,JJ c IF (I.LE.K) GOTO 350 c FO8(1) = FL8(I-K) c L = I c 350 CONTINUE c MO = ((I-1)* (MN2-I))/2 + L c MP = MO + JJ - L c WRITE (IWRITE,FO8) (HNP1(J),J=MO,MP) c 360 CONTINUE c WRITE (IWRITE,3040) c IF (M.LT.MNP1) GOTO 340 C c 370 CONTINUE WRITE (IWRITE,3100) c IF (ABS(IPOLPH).EQ.2) GOTO 380 c **** parallel **** 375 IF (MORE.LT.2 .AND. INAST.NE.0)THEN c REWIND (ITAPE) c REWIND (77) !parallel done earlier ENDIF GOTO 390 C 380 CONTINUE IF (MSKIP.GT.MZSLP) CALL RECOV2('TAPERD',' MZSLP',MZSLP,MSKIP) C 390 CONTINUE IF (MSKIP.LT.IPWINIT) + WRITE (IWRITE,3105) MSKIP,IPWINIT 400 continue !! parallel c if (istop.ne.0) return C 3000 FORMAT (/35X,'SUBROUTINE TAPERD'/35X,17 ('-')) 3020 FORMAT (' CANNOT FIND ON INPUT FILE ANY DATA FOR THIS CASE'/) 3030 FORMAT (' TARGET OR CORE STATE DATA'/8X,' L', A 'PARITY CHANNELS ENERGY(AU) CORRECTION ENERGY(RY)') 3040 FORMAT (A8,I3,1X,A4,I7,F16.7,F12.7,F12.5,I8) 3050 FORMAT (' MAXNHF =',20I3) 3060 FORMAT (' MAXNLG =',20I3) 3070 FORMAT (' NELC =',I3,' NZ =',I3,' LRANG1 =',I3,' LRANG2 =',I3, A ' NRANG2 =',I3,' LAMAX =',I2,' LAM=', B I1/' MASS-CORRECTION(',I1,'), DARWIN-TERM(',I1, C '), SPIN-ORBIT(',I2,');',7X,'IZESP =',I3) 3090 FORMAT (/' HAMILTONIAN MATRIX - SIZE =',I5,I10,' CHANNELS'/) 3100 FORMAT (/' CORRECT DATA READ FROM STG2') 3105 FORMAT (/' (CALCULATION SKIPPED - SYMMETRY=: ',I3,' .LT. ',I3) 3110 FORMAT (' RA =',F10.5,', BSTO =',E12.5,'; NIX =',I2/) 3130 FORMAT (' LRGL =',I3,' NSPN =',I3,' NPTY =',I3,' NCFGP =',I5) 3140 FORMAT ( A ' SPECIFIED CASE OUT OF STEP WITH THIS CASE WHILE IPOLPH.GT.1' B /' REMAINING CASES SKIPPED') 3200 FORMAT('THE FILE ',A12,' DOES NOT EXIST .... ') 3210 FORMAT(' OPENING FILE ',A12/) 3205 FORMAT(' ') 3220 FORMAT('WILL DIAGONALISE PRESENT SYMMETRY, THEN STOP') return END C********************************************************************** C********************************************************************** SUBROUTINE TAPERD(HNP1,irowsize,icolsize) C----------------------------------------------------------------------- C C C C----------------------------------------------------------------------- C C READS THE INPUT DATA FROM STG2 CORRESPONDING TO THE REQUIRED C VALUES OF LRGL, NSPN, NPTY, NCFGP. C C----------------------------------------------------------------------- C c use param use basic2 use boundH use comm_interface IMPLICIT REAL*8 (A-H,O-Z) INTEGER*8 JBIG,JJ,IJJ,IPOS,KP,KQ,KN,MO,MP,LBIG,NQ,MN2BIG,MNP1BIG INTEGER*8 :: L2PBIG(MZCHF) C PARAMETER (MXN21=MZNR2+1) PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) c PARAMETER (MXTRI=(MXMAT*(MXMAT+1))/2) c PARAMETER (MXD1=MZNC2/MZNR2,MXD2=MZNR2/MZNC2,MXD3=MXD1+MXD2, c A MXDM=MZNC2*MXD1/MXD3+MZNR2*MXD2/MXD3) c PARAMETER (MXDMAT=4*MXMAT*MZCHF+2*MXMAT*MXDM+2*MZNR2*MXDM, c A MXRSCT=MXTRI+MZCHF*MXMAT+MZNR2*MXMAT) c PARAMETER (MXM1=MXDMAT/MXRSCT,MXM2=MXRSCT/MXDMAT,MXM3=MXM1+MXM2, c A MXMATX=MXDMAT*MXM1/MXM3+MXRSCT*MXM2/MXM3+1) c PARAMETER (MXDUM2=MXMATX-MXTRI-MZNC2*MZNR2) C LOGICAL EX,EXA,EXBP,EXH CHARACTER PARITY(0:1)*4,SPIN(0:9)*8 CHARACTER*6 FL8(8),FO8(2) CHARACTER*1 IFNM(10) CHARACTER*12 IFLNM CHARACTER*9 IFLBP C COMMON /ADDE/EST(MZTAR),EPOLE,IPOLE,NEST,AST(MZTAR),MERGE(MZTAR) COMMON /BASIC1/BSTO,RA,NELC,NZ,LRANG2,NRANG2,MAXNHF(MZLR2), A MAXNLG(MZLR2),MAXNC(MZLR2),LRANG1 COMMON /BASIN/EIGENS(MZNR2,MZLR2),ENDS(MXN21,MZLR2),DELTA,ETA COMMON /BUTT/COEFF(3,MZLR2),EK2MAX,EK2MIN,MAXNCB(MZLR2),NELCOR COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST,IONEONE COMMON /CUPINT/MNP1,NCONHP,NCHAN,N2HDAT,NTARSYM COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /INFORM/IREAD,IWRITE,IPUNCH c **** parallel **** c COMMON /MATRIX/HNP1(MXTRI),STORE(MZNC2,MZNR2),DUM2(MXDUM2) c **** parallel **** COMMON /MORDER/LAMAX,LAMBC,LAMCC,LAMIND,LAM COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /RCASES/NOT1,NOT2,NOT3,NPROG,NCFGP COMMON /RECOV/IPLACE COMMON /STATE/ENAT(MZTAR),LAT(MZTAR),ISAT(MZTAR),LPTY(MZTAR),NAST COMMON /SHIFT/ESHIFT,NSHIFT COMMON /NRBSRT/ECORR(MZTAR),NCORR,ISORT COMMON /NRBCOD/ICODE COMMON /NRBMNP/MNP2MX,IPRCENT,iprint COMMON /PART/TRACE COMMON /NRBHDX/EXH C DIMENSION JRELOP(3),NORDER(MZTAR) C SAVE IFL,nsym C DATA IFL/1/,nsym/0/ DATA PARITY/'EVEN',' ODD'/,TWO/2.0D0/, A ZERO/0.0D0/,SPIN/'(= 2J) ',' SINGLET', A ' DOUBLET',' TRIPLET',' QUARTET',' QUINTET',' SEXTET ', B ' SEPTET ',' OCTET ',' 2S+1>8 '/ DATA FL8/' (8','(14X,7','(28X,6','(42X,5','(56X,4','(70X,3', A '(84X,2',' (98X,'/,FO8(2)/'F14.7)'/ DATA IFNM/'0','1','2','3','4','5','6','7','8','9'/ c **** parallel **** common /crashblock/istop include 'mpif.h' common /procdat/ictxt,myrow,mycol c common /parablock/iam,nproc dimension HNP1(irowsize,icolsize) allocatable store(:,:) allocatable haux(:) allocatable hcorr(:) COMMON /REDTOL/TOLER common /contsize/icontold,nonzer common /globsize/idimH(mzslp) character*1 filec,filed,fileu common /allocblock/iflagallhbb,iflagallhbc,iflagallcf c **** parallel **** C----------------------------------------------------------------------- C C IF ANY DIMENSION IS EXCEEDED WHEN READING FROM DATA, CALL RECOV2 C WITH IPLACE=0 TO TERMINATE THE PROGRAM C IPLACE = 0 WRITE (IWRITE,3000) ITAPE = ITAPE2 TRACE=ZERO !partitioned MSKIP0=MSKIP c nsym=nsym+1 c write(0,*)'mskip,nsym,ipolh=',mskip,nsym,IPOLPH inquire(file='sizeBP.dat',exist=exbp) if(.not.exbp.or.nsym.gt.1)then IF (MSKIP.GT.1 .AND. (INAST.EQ.0.OR. A ABS(IPOLPH).EQ.2)) GOTO 20 else IF (MSKIP.GT.1 .AND. (INAST.EQ.0.OR. A ABS(IPOLPH).EQ.2)) GOTO 8 endif c **** parallel **** C REWIND (77) !already called c **** parallel **** REWIND (ITAPE) C C READ THE BASIC DATA FROM FILE C 8 IF(N2HDAT.GT.0) THEN READ (ITAPE) NELC,NZ,LRANG1,LRANG2,NRANG2,LAMAX,ICODE,LAM,IZESP, A (JRELOP(I),I=1,3),MAXLD ELSE READ (ITAPE,ERR=9) NELC,NZ,LRANG1,LRANG2,NRANG2,LAMAX,ICODE,LAM, A IZESP,(JRELOP(I),I=1,3),MAXLD,N2HDT0 N2HDAT=N2HDT0 END IF 9 WRITE (IWRITE,3070) NELC,NZ,LRANG1, A LRANG2,NRANG2,LAMAX,LAM,JRELOP,IZESP N2HDAT=-1 IF (LRANG1.GT.MZLR2) CALL RECOV2('TAPERD',' MZLR2',MZLR2,LRANG1) IF (LRANG2.GT.MZLR2) CALL RECOV2('TAPERD',' MZLR2',MZLR2,LRANG2) IF (NRANG2.GT.MZNR2) CALL RECOV2('TAPERD',' MZNR2',MZNR2,NRANG2) IF (LAMAX.GT.MZLMX) CALL RECOV2('TAPERD',' MZLMX',MZLMX,LAMAX) if (istop.ne.0) return NOT2 = MIN(NRANG2,NOT2) !no longer used READ (ITAPE) (MAXNHF(L),L=1,LRANG1), (MAXNLG(L),L=1,LRANG1), A (MAXNC(L),L=1,LRANG1) WRITE (IWRITE,3050) (MAXNHF(L),L=1,LRANG1) WRITE (IWRITE,3060) (MAXNLG(L),L=1,LRANG1) DO 10 L = 1,LRANG2 READ (ITAPE) (EIGENS(N,L),N=1,NRANG2) READ (ITAPE) (ENDS(N,L),N=1,NRANG2) 10 CONTINUE READ (ITAPE) RA,BSTO,HINT,DELTA,ETA,NIX WRITE (IWRITE,3110) RA,BSTO,NIX C C IF(NIX.GT.0) THEN C READ(ITAPE) (IHX,I=1,NIX),(IRX,I=1,NIX) C IPTS=2*IRX C READ(ITAPE) (POVALU,I=1,IPTS) C DO 3 LP=1,LRANG1 C NBT=MAXNLG(LP)-LP+1 C IF(NBT.LE.0) GOTO 2 C DO 2 K=1,NBT C 2 READ(ITAPE) (POVALU,I=1,IPTS) C 3 CONTINUE C ENDIF if(exbp.AND.ifl.ne.1)then !Since RECUPDXXX complete, but don't want READ (ITAPE) !to overwrite ordered/observed energies READ (ITAPE) READ (ITAPE) goto 20 endif READ (ITAPE) ((COEFF(I,L),I=1,3),L=1,LRANG2) READ (ITAPE) NAST WRITE(IWRITE,*)' NAST = ',NAST NAST = ABS(NAST) IF (NAST.GT.MZTAR) CALL RECOV2('TAPERD',' MZTAR',MZTAR,NAST) if (istop.ne.0) return READ (ITAPE) (ENAT(I),I=1,NAST), (LAT(I),I=1,NAST), A (ISAT(I),I=1,NAST), (LPTY(I),I=1,NAST) C C THE FOLLOWING DATA DEPENDS ON LRGL, NSPN ETC. C 20 CONTINUE ieof=1 READ (ITAPE,END=210) LRGL2,NSPN2,NPTY2,NCFGP2,IPOL2 WRITE (IWRITE,3130) LRGL2,NSPN2,NPTY2,NCFGP2 ieof=0 C READ (ITAPE) MNP1,NCONHP,NCHAN c if(iam.eq.0)write(*,*)mnp1,nconhp,nchan,mskip if(idimH(mskip).ne.mnp1)then !file mis-match if(iam.eq.0)then write(iwrite,*)'For symmetry=',mskip,' idimH=',idimH(mskip) x ,' but H-file read mnp1=',mnp1,' for SLp=',nspn2,lrgl2,npty2 if(ioneone.gt.0)write(iwrite,*)'Incorrect ioneone=',ioneone x ,', the number of symmetries in each H-file?' write(0,*)'Mis-match between expected and found' x ,' on file symmetry:',mskip endif call flush(iwrite) stop !need something drastic else job just hangs c istop=42 c return endif C IF(IONEONE.gt.0.AND.MSKIP-1+IONEONE.LT.IPWINIT)THEN WRITE(IWRITE,*)' ' !'MSKIP= ',MSKIP WRITE(IWRITE,*)'SKIPPING EARLY STG2HXXX.DAT IF EARLY SYMMETRIES' MSKIP=MSKIP+IONEONE-1 WRITE(IWRITE,*)'ARE NOT REQUIRED and WE HAVE ' X ,IONEONE,' SYM PER FILE' WRITE(IWRITE,*)' ' !'MSKIP= ',MSKIP c call flush(iwrite) GOTO 210 ENDIF C IF (NCHAN.GT.MZCHF) CALL RECOV2('TAPERD',' MZCHF',MZCHF,NCHAN) IF (MNP1.GT.MXMAT) CALL RECOV2('TAPERD','MXMAT ',MXMAT,MNP1) IF (NCFGP2.GT.MZNC2) CALL RECOV2('TAPERD',' MZNC2',MZNC2,NCFGP2) if (istop.ne.0) return READ (ITAPE) (NCONAT(I),I=1,NAST) IF(NSPN2.NE.0)READ (ITAPE) (L2P(I),I=1,NCHAN) IF(NSPN2.EQ.0)READ (ITAPE) (L2P(I),I=1,NCHAN),(KJ(I),I=1,NCHAN) ! NRB READ (ITAPE) MORE2 c if(iam.eq.0)write(*,*)more2 c call flush_(iwrite) c c if HDCORR.DAT exists, subtract corrections from diagonal c IF(EXH.AND.NCFGP2.GT.0)THEN READ(31,*)LRGL,NSPN,NPTY,mdiag IF(LRGL.NE.LRGL2.OR.NSPN.NE.NSPN2.OR.NPTY.NE.NPTY2)THEN WRITE(IWRITE,*) 'MIS-MATCH ON HDCORR.DAT FOR L S P:', X LRGL2,NSPN2,NPTY2,LRGL,NSPN,NPTY STOP 'MIS-MATCH ON HDCORR.DAT FOR L S P' ENDIF IF(MNP1.NE.mdiag)THEN WRITE(IWRITE,*) 'MIS-MATCH ON HDCORR.DAT FOR MNP1:', X MNP1,mdiag STOP 'MIS-MATCH ON HDCORR.DAT FOR MNP1' ENDIF c **** parallel **** allocate(hcorr(mnp1),stat=ierr) if(ierr.ne.0)stop'program could not allocate memory for hcorr' c **** parallel **** do j=1,MNP1 read(31,*)jjj,hcorr(j) enddo ENDIF C c **** parallel **** allocate(store(nrang2,nrang2),stat=ierr) if(ierr.ne.0)stop'program could not allocate memory for store' c **** parallel **** C C READ THE UPPER TRIANGLE OF THE HAMILTONIAN MATRIX INTO HNP1 C do M=1,NCHAN L2PBIG(M)=L2P(M) enddo C MDIAGCOUNT=1 !partitioned ICOUNTA=0 !partitioned mdiag=0 !hcorr C c if(iam.eq.0)write(*,*)'starting CC read',icode MN2 = 2*MNP1 MN2BIG=MN2 MNP1BIG=MNP1 C DO 110 M = 1,NCHAN MP = (M* (M+1))/2 MO = MP - M + 1 NQ = (M-1)*NRANG2 NBT = 1 DO 100 LBIG = MO,MP C C ---> READ CONTINUUM-CONTINUUM BLOCKS INTO STORE IN TURN IF (ICODE.NE.21) GOTO 50 IF (LBIG.EQ.MP) GOTO 40 IF (L2PBIG(M).EQ.L2PBIG(LBIG-MO+1)) GOTO 50 40 CONTINUE if(iam.eq.0.and.MSKIP.GE.IPWINIT)then READ (ITAPE) ((STORE(I,J),I=1,NRANG2),J=1,NRANG2) else READ (ITAPE) endif IF (LBIG.GE.MP) GOTO 70 GOTO 60 C 50 CONTINUE if(iam.eq.0.and.MSKIP.GE.IPWINIT)then READ (ITAPE) ((STORE(J,I),I=1,NRANG2),J=1,NRANG2) else READ (ITAPE) endif 60 CONTINUE KP = (LBIG-MO)*NRANG2 KN = ((MN2BIG-KP+1)*KP)/2 + NQ - KP C C C ---> TRANSFER CONTINUUM-CONTINUUM BLOCKS INTO HNP1 FROM STORE 70 CONTINUE C C SKIP IF THIS SYMMETRY IS NOT NEEDED IF (MSKIP.LT.IPWINIT) GO TO 100 c n222=nrang2*nrang2 call mpi_bcast(store,n222,mpi_real8,0,new_comm,ierr) C DO 90 J = 1,NRANG2 IF (LBIG.NE.MP) THEN JBIG=J JJ = (JBIG-1)* (MNP1BIG-KP) - ((JBIG-1)*JBIG)/2 + KN C ELSE NBT = J JBIG=J JJ = ((NQ+JBIG-1)* (MN2BIG-NQ-JBIG+2))/2 - JBIG + 1 ENDIF C DO 80 I = NBT,NRANG2 c **** parallel **** c...........in serial stg3: HNP1(I+JJ) = STORE(I,J) c...........here I+JJ global position is mapping into local (irow,icol) if (toler.lt.1.d-19.or.nonzer.eq.0) then ipos=I+JJ else c...............transform the index I+JJ for size MNP1 to c...............the index ipos for size idimH(mskip) IJJ=I+JJ c call ntoij(mnp1,IJJ,ig,jg) c call ijton(idimH(mskip),ipos,ig,jg) call ntoijandback(mnp1,IJJ,ipos) endif C C Locate diagonal both for HCC contribution to TRACE & hcorr C if(IPRCENT.ne.100.or.EXH)then if(ipos.eq.mdiagcount)then if(EXH.AND.NCFGP2.GT.0)then mdiag=mdiag+1 store(i,j)=store(i,j)-hcorr(mdiag) endif if(IPRCENT.ne.100)TRACE=TRACE+store(i,j) !Partitioned mdiagcount=mdiagcount+(mnp1-icounta) icounta=icounta+1 endif endif c call locH(ipos,irow,icol,ipr,ipc) if (ipr.eq.myrow.and.ipc.eq.mycol) + HNP1(irow,icol) = store(i,j) C c **** parallel **** 80 CONTINUE 90 CONTINUE 100 CONTINUE 110 CONTINUE c 111 continue c deallocate(store) C c if(iam.eq.0)write(*,*)'ending CC read' C C READ BOUND-CONTINUUM BLOCKS INTO STORE C IF (NCFGP2.EQ.0) GOTO 160 IF (NCFGP2.GT.MZNC2) CALL RECOV2('TAPERD',' MZNC2',MZNC2,NCFGP2) if (istop.ne.0) return c **** parallel **** if (toler.ge.1.d-19.and.nonzer.gt.0.or.nshift.gt.0) then if (iflagallhbc.ne.0) then deallocate(hbc) iflagallhbc = 0 endif if (iflagallhbb.ne.0) then deallocate(hbb) iflagallhbb = 0 endif if(toler.ge.1.d-19.and.nonzer.gt.0)then ncont=mnp1-ncfgp2 allocate(hbc(ncfgp2,ncont),stat=ierr) if(ierr.ne.0)stop'allocation fails for hbc' iflagallhbc = 1 endif allocate(hbb(ncfgp2,ncfgp2),stat=ierr) if(ierr.ne.0)stop'allocation fails for hbb' iflagallhbb = 1 endif allocate(store(ncfgp2,nrang2),stat=ierr) if(ierr.ne.0) stop 'allocation fails for store (b)' c **** parallel **** c if(iam.eq.0)write(*,*)'starting BC read' C DO 140 M = 1,NCHAN C C SKIP IF THIS SYMMETRY IS NOT NEEDED IF (MSKIP.LT.IPWINIT) THEN READ (ITAPE) GO TO 140 ENDIF C READ (ITAPE) ((STORE(I,J),I=1,NCFGP2),J=1,NRANG2) C KP = (M-1)*NRANG2 KQ = (KP* (MN2BIG-KP+1))/2 + NCONHP DO 130 J = 1,NRANG2 JBIG=J JJ = (JBIG-1)*MNP1BIG - (JBIG* (JBIG-1))/2 - JBIG*KP + KQ DO 120 I = 1,NCFGP2 c **** parallel **** c........in serial stg3: HNP1(I+JJ) = STORE(I,J) if (toler.lt.1.d-19.or.nonzer.eq.0) then c............here I+JJ global position is mapping into local (irow,icol) IJJ=I+JJ call locH(IJJ,irow,icol,ipr,ipc) if (ipr.eq.myrow.and.ipc.eq.mycol) + HNP1(irow,icol) = store(i,j) else c...............put in HBC for further transformation in HAMNEW IJJ=I+JJ call ntoij(mnp1,IJJ,irow,icol) icol=icol-icontold c................HBC is the transverse of HCB HBC(icol,irow) = store(i,j) endif c **** parallel **** 120 CONTINUE 130 CONTINUE 140 CONTINUE c 141 continue c deallocate(store) allocate(haux(ncfgp2),stat=ierr) if(ierr.ne.0)stop'allocation fails for haux' C c if(iam.eq.0)write(*,*)'ending BC read' C C READ BOUND-BOUND ELEMENTS INTO HNP1 C c if(iam.eq.0)write(*,*)'starting BB read' c if(mdiag.ne.nconhp)stop 'mdiag.ne.nconhp' C MP = (NCONHP* (MN2BIG-NCONHP+1))/2 C c if(npw_per_subworld.eq.1.and.mskip.ne.(icolour+1))then c goto 151 c endif C DO 150 I = 1,NCFGP2 MO = MP + 1 MP = MO + NCFGP2 - I c **** parallel **** c........in serial stg3: READ (ITAPE) (HNP1(J),J=MO,MP) c........here an auxiliar haux array is created which has hnp1(J) J=mo,mp c........then, the J global position is mapping into local (irow,icol) c c........skip if this symmetry is not needed if (mskip.lt.ipwinit) then read (itape) go to 150 endif c read (itape) (haux(j),j=1,mp-mo+1) c jj=1 C Locate diagonal both for HBB contribution to TRACE & hcorr if(EXH.AND.NCFGP2.GT.0)haux(jj)=haux(jj)-hcorr(nconhp+i) if(IPRCENT.ne.100)TRACE=TRACE+haux(jj) !Partitioned c C do 149 j=mo,mp do 149 jbig=mo,mp c if ((toler.lt.1.d-19.or.nonzer.eq.0).and.nshift.le.0) then c call locH(j,irow,icol,ipr,ipc) call locH(JBIG,irow,icol,ipr,ipc) if (ipr.eq.myrow.and.ipc.eq.mycol) + HNP1(irow,icol) = haux(jj) else c.............put in HBB for further transformation in HAMNEW(2) c call ntoij(mnp1,J,irow,icol) call ntoij(mnp1,jbig,irow,icol) irow=irow-icontold icol=icol-icontold HBB(irow,icol) = haux(jj) HBB(icol,irow) = haux(jj) !symmetrize endif jj=jj+1 149 continue c **** parallel **** 150 CONTINUE c 151 continue C c deallocate(haux) if(EXH.AND.NCFGP2.GT.0)deallocate(hcorr) C c if(iam.eq.0)write(*,*)'ending BB read' C C READ THE ASYMPTOTIC COEFFICIENTS (ARRAY CF EXPANDED IN WRITOP) C 160 CONTINUE IF (LAMAX.EQ.0) GOTO 180 c c **** parallel **** if (mskip.ge.ipwinit) then if (iflagallcf.ne.0) then !TBD skip deallocate(cf) iflagallcf = 0 endif allocate(cf(nchan,nchan,lamax),stat=ierr) if(ierr.ne.0)stop'allocation fails for cf' iflagallcf = 1 endif c **** parallel **** c IF (MSKIP.LT.IPWINIT) THEN DO 170 I = 1,NCHAN READ (ITAPE) 170 CONTINUE ELSE DO 175 I = 1,NCHAN READ (ITAPE) ((CF(I,J,K),K=1,LAMAX),J=I,NCHAN) 175 CONTINUE ENDIF C C CHECK THAT LRGL,NSPN,NPTY,NCFGP HAVE THE REQUIRED VALUES. C IF NOT, READ IN MORE DATA FROM FILE OR USER. C 180 CONTINUE IF (MORE2.EQ.0.and..not.exbp) THEN MORE = 0 END IF IF (INAST.NE.0) GOTO 190 LRGL = LRGL2 NSPN = NSPN2 NPTY = NPTY2 190 CONTINUE c IF (LRGL.EQ.LRGL2 .AND. NSPN.EQ.NSPN2 .AND. A NPTY.EQ.NPTY2) THEN if(exbp.and.more2.eq.0)goto 210 GOTO 230 ENDIF c IF (IPOLPH.LE.1) GOTO 200 WRITE (IWRITE,3140) MORE = 0 GOTO 220 C 200 CONTINUE IF (MORE2.GE.1) THEN if (iflagallcf.ne.0) then deallocate(cf) iflagallcf = 0 endif if(exbp)go to 8 GOTO 20 ENDIF 210 CONTINUE C IF(N2HDAT.EQ.IFL) THEN !NO MORE FILES IF(INAST.EQ.0) THEN MORE=0 ELSE WRITE (IWRITE,3020) IF(MSKIP.LT.INAST) MORE=1 END IF ELSE C C OPEN NEW FILE STG2HXXX.DAT ich0 = ichar('0') ic = ifl/100 id = (ifl - 100*ic)/10 iu = (ifl - 100*ic - 10*id) ic = ic + ich0 id = id + ich0 iu = iu + ich0 filec = char(ic) filed = char(id) fileu = char(iu) c nsym=0 c if(.not.exbp)then c IFLNM='STG2H'//filec//filed//fileu//'.DAT' INQUIRE (FILE=IFLNM,EXIST=EX) IF(.NOT.EX) THEN WRITE(IWRITE,3200) IFLNM istop=4 return ELSE CLOSE(ITAPE) OPEN (UNIT=ITAPE,FILE=IFLNM,STATUS='OLD',FORM='UNFORMATTED') WRITE(IWRITE,3210) IFLNM END IF IFLNM='HDCORR'//filec//filed//fileu INQUIRE (FILE=IFLNM,EXIST=EXH) IF(EXH)THEN CLOSE(31) OPEN (31,FILE=IFLNM,STATUS='OLD',FORM='FORMATTED') WRITE(IWRITE,3210) IFLNM ENDIF IFL=IFL+1 if(ieof.eq.1.and.mskip.gt.1)GO TO 20 GO TO 230 c else c IFLBP='RECUPH'//filec//filed//fileu INQUIRE (FILE=IFLBP,EXIST=EX) IF(.NOT.EX) THEN WRITE(IWRITE,3205) WRITE(IWRITE,3200) IFLBP WRITE(IWRITE,3205) WRITE(IWRITE,3220) WRITE(IWRITE,3205) MORE=0 ELSE CLOSE(ITAPE) OPEN (UNIT=ITAPE,FILE=IFLBP,STATUS='OLD',FORM='UNFORMATTED') WRITE(IWRITE,3205) WRITE(IWRITE,3210) IFLBP WRITE(IWRITE,3205) MORE=1 END IF IFL=IFL+1 c IPLACE = -1 GOTO 230 c endif c ENDIF c 220 CONTINUE IPLACE = -1 GOTO 390 230 CONTINUE C C SKIP IF THIS SYMMETRY IS NOT NEEDED IF (MSKIP.LT.IPWINIT) GO TO 235 C C C TRANSFORM THE HAMILTONIAN TO REDUCE NUMBER OF N+1-CFGS - TWG 13/4/95. C IF(NCFGP2.GT.0.and.TOLER.ge.1.D-19.and.nonzer.gt.0)THEN call HAMNEW(HNP1,irowsize,icolsize, + MNP1,NCFGP2,LRGL2,NSPN2,NPTY2) MN2=2*MNP1 ENDIF NCFGP = NCFGP2 C C SHIFT THE N+1 CFGS - TWG 18/11/95. C IF(NCFGP2.GT.0.and.NSHIFT.gt.0) THEN !parallel CALL HAMNEW2(HNP1,irowsize,icolsize, + MNP1,NCFGP2,LRGL2,NSPN2,NPTY2) ENDIF c NCFGP = NCFGP2 !not redefined c **** parallel **** 235 if(iflagallhbc.ne.0)then deallocate(hbc) iflagallhbc = 0 endif if(iflagallhbb.ne.0)then deallocate(hbb) iflagallhbb = 0 endif c **** parallel **** C C WRITE OUT THE TARGET STATES. C IF (MSKIP0.GT.1) GOTO 250 WRITE (IWRITE,3030) CALL ORDER(ENAT,NORDER,NAST,1) IF (NCORR+ABS(NEST).EQ.0) THEN GROUND = ENAT(NORDER(1)) DO 240 II = 1,NAST I = NORDER(II) T = (ENAT(I)-GROUND)*2 IS = MIN(ISAT(I),9) WRITE (IWRITE,3040) SPIN(IS),LAT(I), A PARITY(LPTY(I)),NCONAT(I),ENAT(I),ZERO,T,II 240 CONTINUE ENDIF C C IF NEST=NAST, ADD ON THE EXTRA ENERGY TO EACH STATE AND TO THE C ASSOCIATED DIAGONAL ELEMENTS OF THE HAMILTONIAN MATRIX. C C EXTENSION '90JUL1/5: OBS.ENERGY INPUT IN RYD ABOVE 1 OR IN CM**-1 C EST(I)=0 ENSURES UNCORRECTED TERM ENERGY IF NO OBS.DATA AVAILABLE C 250 CONTINUE IF (NCORR+ABS(NEST).EQ.0) GOTO 310 E0 = ENAT(1) IF (MSKIP0.LE.1) THEN C C IF NEST .LT. 0 GROUP ENERGIES AROUND ARITHMETIC MEANS C IF (NEST.LT.0) THEN DO II = 1,NAST EST(II) = ENAT(II) ENDDO ISTART = 0 DO 258 I = 1,-NEST IF(ISTART.GE.NAST.OR.MERGE(I).LE.0) GOTO 258 IF(MERGE(I).EQ.1) GOTO 256 EMEAN = 0.0D0 IEND = MIN(ISTART+MERGE(I),NAST) DO II=ISTART+1,IEND EMEAN = EMEAN+ENAT(II) ENDDO EMEAN = EMEAN / (IEND-ISTART) DO II=ISTART+1,IEND EST(II) = EMEAN ENDDO 256 ISTART = ISTART + MERGE(I) 258 CONTINUE C IF (ISTART.NE.NAST) THEN WRITE(IWRITE,*) '**** WARNING **** ', A 'Target State Energy MERGE only treats', B ISTART,' states' ENDIF ENDIF C C IF NCORR,ABS(NEST) .GT. 0 CORRECT ENERGIES C IF(EST(1).EQ.ZERO)E0=E0+ESHIFT DO 270 II = 1,NAST I = NORDER(II) I0=I IF(ISORT.GT.0)I0=II C AST(I) = ZERO IF(NEST.NE.0)AST(I) = EST(I0) - ENAT(I) IF(NCORR.GT.0)AST(I)=AST(I)+ECORR(I0) C IF (EST(1).EQ.ZERO) THEN AST(I) = ZERO IF(NEST.NE.0)AST(I) = EST(I0) IF(NCORR.GT.0)THEN IF(ECORR(I0).NE.0)THEN AST(I)=AST(I)+ECORR(I0) IF(EST(I0).EQ.0)AST(I)=AST(I)+(ENAT(I)-E0)*TWO ENDIF ENDIF AST(I)=AST(I)/TWO C IF (AST(I).EQ.ZERO)THEN ENAT(I) = ENAT(I) + ESHIFT GOTO 260 ENDIF C IF (EST(I0).LT.ZERO) AST(I) = -AST(I)/109737.3D0 AST(I) = AST(I) - ENAT(I) + E0 ENDIF C ENAT(I) = ENAT(I) + AST(I) C 260 CONTINUE T = (ENAT(I)-ENAT(1))*TWO IS = MIN(ISAT(I),9) WRITE(IWRITE,3040)SPIN(IS),LAT(I),PARITY(LPTY(I)),NCONAT(I), A ENAT(I),AST(I),T,II 270 CONTINUE ELSE IF(INAST.GT.0)THEN C ENAT HAS BEEN RE-READ AND SO MUST BE RE-ADJUSTED : NRB 4/11/94 CALL ORDER(ENAT,NORDER,NAST,1) DO 275 II=1,NAST I=NORDER(II) ENAT(I)=ENAT(I)+AST(I) 275 CONTINUE ENDIF C END MODS ENDIF C C C LOOP OVER THE CHANNELS ASSOCIATED WITH EACH STATE IN THE H-MATRIX C C SKIP IF THIS SYMMETRY IS NOT NEEDED IF (NEST.EQ.0.OR.MSKIP.LT.IPWINIT) GOTO 310 C MN2BIG=MN2 !CPB oct 2005 (as MN2 has been redefined) C MP = 0 DO 300 I = 1,NAST IF (NCONAT(I).LE.0) GOTO 300 MO = MP + 1 MP = MP + NCONAT(I) DO 290 M = MO,MP NQ = (M-1)*NRANG2 DO 280 J = 1,NRANG2 JBIG=J JJ = ((NQ+JBIG-1)* (MN2BIG-NQ-JBIG+2))/2 + 1 c **** parallel **** c............in serial stg3: HNP1(JJ) = HNP1(JJ) + AST(I) call locH(JJ,irow,icol,ipr,ipc) if (ipr.eq.myrow.and.ipc.eq.mycol) + HNP1(irow,icol) = HNP1(irow,icol) + AST(I) c **** parallel **** 280 CONTINUE 290 CONTINUE 300 CONTINUE C C WRITE THE FULL HAMILTONIAN MATRIX TO SCRATCH DISC C 310 CONTINUE C SKIP IF THIS SYMMETRY WAS NOT READ IF (IONEONE.GT.0.AND.MSKIP.LT.IPWINIT) GO TO 390 c **** parallel **** c all this chunk is commented - need to be implemented in the future c IF (NOT1.GE.NRANG2) GOTO 330 c REWIND IDISC1 c DO 320 M = 1,NCHAN c MO = ((M-1)*NRANG2* (MN2- (M-1)*NRANG2+1))/2 + 1 c MP = (M*NRANG2* (MN2-M*NRANG2+1))/2 c WRITE (IDISC1) (HNP1(I),I=MO,MP) c 320 CONTINUE c MO = (NCONHP* (MN2-NCONHP+1))/2 + 1 c MP = (MNP1* (MNP1+1))/2 c WRITE (IDISC1) (HNP1(I),I=MO,MP) C C WRITE OUT THE HAMILTONIAN MATRIX IF NBUG6=1 C c 330 CONTINUE WRITE (IWRITE,3090) MNP1,NCHAN c IF (NBUG6.NE.1) GOTO 370 c M = 0 c 340 CONTINUE c K = M c L = M + 1 c M = M + 8 c JJ = MIN(M,MNP1) c FO8(1) = FL8(1) c DO 360 I = 1,JJ c IF (I.LE.K) GOTO 350 c FO8(1) = FL8(I-K) c L = I c 350 CONTINUE c MO = ((I-1)* (MN2-I))/2 + L c MP = MO + JJ - L c WRITE (IWRITE,FO8) (HNP1(J),J=MO,MP) c 360 CONTINUE c WRITE (IWRITE,3040) c IF (M.LT.MNP1) GOTO 340 C c 370 CONTINUE WRITE (IWRITE,3100) c IF (ABS(IPOLPH).EQ.2) GOTO 380 c **** parallel **** 375 IF (MORE.LT.2 .AND. INAST.NE.0)THEN REWIND (ITAPE) c REWIND (77) !parallel done earlier ENDIF GOTO 390 C 380 CONTINUE IF (MSKIP.GT.MZSLP) CALL RECOV2('TAPERD',' MZSLP',MZSLP,MSKIP) C 390 CONTINUE IF (MSKIP.LT.IPWINIT) + WRITE (IWRITE,3105) MSKIP,IPWINIT 400 continue !! parallel c if (istop.ne.0) return C 3000 FORMAT (/35X,'SUBROUTINE TAPERD'/35X,17 ('-')) 3020 FORMAT (' CANNOT FIND ON INPUT FILE ANY DATA FOR THIS CASE'/) 3030 FORMAT (' TARGET OR CORE STATE DATA'/8X,' L', A 'PARITY CHANNELS ENERGY(AU) CORRECTION ENERGY(RY)') 3040 FORMAT (A8,I3,1X,A4,I7,F16.7,F12.7,F12.5,I8) 3050 FORMAT (' MAXNHF =',20I3) 3060 FORMAT (' MAXNLG =',20I3) 3070 FORMAT (' NELC =',I3,' NZ =',I3,' LRANG1 =',I3,' LRANG2 =',I3, A ' NRANG2 =',I3,' LAMAX =',I2,' LAM=', B I1/' MASS-CORRECTION(',I1,'), DARWIN-TERM(',I1, C '), SPIN-ORBIT(',I2,');',7X,'IZESP =',I3) 3090 FORMAT (/' HAMILTONIAN MATRIX - SIZE =',I5,I10,' CHANNELS'/) 3100 FORMAT (/' CORRECT DATA READ FROM STG2') 3105 FORMAT (/' (CALCULATION SKIPPED - SYMMETRY=: ',I3,' .LT. ',I3) 3110 FORMAT (' RA =',F10.5,', BSTO =',E12.5,'; NIX =',I2/) 3130 FORMAT (' LRGL =',I3,' NSPN =',I3,' NPTY =',I3,' NCFGP =',I5) 3140 FORMAT ( A ' SPECIFIED CASE OUT OF STEP WITH THIS CASE WHILE IPOLPH.GT.1' B /' REMAINING CASES SKIPPED') 3200 FORMAT('THE FILE ',A12,' DOES NOT EXIST .... ') 3210 FORMAT(' OPENING FILE ',A12/) 3205 FORMAT(' ') 3220 FORMAT('WILL DIAGONALISE PRESENT SYMMETRY, THEN STOP') return END C********************************************************************** SUBROUTINE VMUL(A,B,C,L,M,N) c use param IMPLICIT REAL*8 (A-H,O-Z) C C C C----------------------------------------------------------------------- C C MATRIX MULTIPLICATION A*B = C, C OUTPUTS RESULTS FOR C(1:L,1:N) C C----------------------------------------------------------------------- C PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) C DIMENSION A(MZNR2,MXMAT),B(MXMAT,N),C(MZNR2,N) C----------------------------------------------------------------------- DO 40 J = 1,N DO 10 I = 1,L C(I,J) = 0.0D0 10 CONTINUE DO 30 K = 1,M DO 20 I = 1,L C(I,J) = A(I,K)*B(K,J) + C(I,J) 20 CONTINUE 30 CONTINUE 40 CONTINUE C END C********************************************************************** SUBROUTINE WRITOP C C C C----------------------------------------------------------------------- C C WRITES THE BASIC QUANTITIES TO THE STG3 FILE C C----------------------------------------------------------------------- C use param use basic2 use comm_interface IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (MXN21=MZNR2+1) C DIMENSION ETEM(MZTAR),NTEM(MZTAR),KTEM(MZTAR),NORDER(MZTAR), A IORDER(MZCHF) C COMMON /BASIC1/BSTO,RA,NELC,NZ,LRANG2,NRANG2,MAXNHF(MZLR2), A MAXNLG(MZLR2),MAXNC(MZLR2),LRANG1 COMMON /BASIN/EIGENS(MZNR2,MZLR2),ENDS(MXN21,MZLR2),DELTA,ETA COMMON /BUTT/COEFF(3,MZLR2),EK2MAX,EK2MIN,MAXNCB(MZLR2),NELCOR COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST,IONEONE COMMON /CUPINT/MNP1,NCONHP,NCHAN,N2HDAT,NTARSYM COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /MORDER/LAMAX,LAMBC,LAMCC,LAMIND,LAM COMMON /RCASES/NOT1,NOT2,NOT3,NPROG,NCFGP COMMON /STATE/ENAT(MZTAR),LAT(MZTAR),ISAT(MZTAR),LPTY(MZTAR),NAST COMMON /NRBCOD/ICODE C COMMON /NRBDIP/MAXLD COMMON /NRBMNP/MNP2MX,IPRCENT,iprint c **** parallel **** include 'mpif.h' c common /parablock/iam,nproc common /matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest X ,NDIV common /allocblock/iflagallhbb,iflagallhbc,iflagallcf C----------------------------------------------------------------------- C c if (lamax.gt.0) then if (iam.ne.0.or.nproctest.ne.0) then if (iflagallcf.ne.0) then deallocate(cf) iflagallcf = 0 endif return endif c endif c **** parallel **** C CALL ORDER(ENAT,NORDER,NAST,1) C ITAPE = ITAPE3 NOTERM=MIN(NOT1,NRANG2) C C FULL H.DAT writes for icomplete_hdat.eq.1 CPB dec 2009 C if(icomplete_hdat.eq.1)goto 5 C C IF MSKIP.GT.1, STG3 IS BEING REPEATED WITH NEW LRGL, NSPN, ETC. C SO IT IS UNNECESSARY TO WRITE OUT THIS FIRST SET OF DATA AGAIN. C IF (MSKIP.GT.1) GOTO 20 C 5 CONTINUE C C REWIND ITAPE LRNG2=LRANG2 !NRB IF(ICODE.EQ.25)LRNG2=-LRANG2 !FLAG FROM DARC NASTY=NAST IF(IPRCENT.NE.100)NASTY=-NAST !FLAG PARTITIONED WRITE (ITAPE) NELC,NZ,LRNG2,LAMAX,NASTY,RA,BSTO C ,MAXLD DO 10 I = 1,NAST I1 = NORDER(I) ETEM(I) = ENAT(I1) KTEM(I) = LAT(I1) NTEM(I) = ISAT(I1) 10 CONTINUE WRITE (ITAPE) (ETEM(I),I=1,NAST) WRITE (ITAPE) (KTEM(I),I=1,NAST) WRITE (ITAPE) (NTEM(I),I=1,NAST) WRITE (ITAPE) ((COEFF(I,L),I=1,3),L=1,LRANG2) C IF(IPRCENT.NE.100)THEN !FOR PARTITIONED WRITE(ITAPE)IPRCENT,NOTERM DO L = 1,LRANG2 WRITE (ITAPE) (EIGENS(N,L),N=1,NOTERM) WRITE (ITAPE) (ENDS(N,L),N=1,NOTERM) ENDDO ENDIF C C THE FOLLOWING DATA DEPENDS ON LRGL, NSPN, NPTY C 20 CONTINUE ETOT = 0.0D0 NF = 0 DO 40 I = 1,NAST IF (NCONAT(I).LE.0) GOTO 40 NS = NF + 1 NF = NF + NCONAT(I) DO 30 K = NS,NF ET(K) = (ETOT-ENAT(I)) !A.U. C *2.0D0 !RYD C ENAT IS A.U. SO IF ET IS RYD THEN THERE CAN BE AN INCONSISTANCY BETWEEN C TERM AND CHANNEL ORDERING BECAUSE OF THE 1.E-7 TOLERANCE IN SR.ORDER - NRB C 30 CONTINUE 40 CONTINUE C C MNP2 = NOTERM*NCHAN + NCFGP MNP2N = (NOTERM*IPRCENT)/100 MNP2N = MNP2N*NCHAN + NCFGP C icomplete_hdat=0 !why do we need this??? C C if(NDIV.eq.1)then WRITE (ITAPE) LRGL,NSPN,NPTY,NCHAN,MNP2N,MORE else WRITE (ITAPE) LRGL,NSPN,NPTY,NCHAN,-MNP2N,MORE endif C CALL ORDER(ET,IORDER,NCHAN,-1) WRITE (ITAPE) (NCONAT(NORDER(I)),I=1,NAST) IF(NSPN.NE.0)WRITE (ITAPE) (L2P(IORDER(I)),I=1,NCHAN) IF(NSPN.EQ.0)WRITE (ITAPE) (L2P(IORDER(I)),I=1,NCHAN) X ,(KJ(IORDER(I)),I=1,NCHAN) ! NRB IF (LAMAX.EQ.0) RETURN DO 70 K = 1,LAMAX DO 60 I = 1,NCHAN DO 50 J = 1,NCHAN CF(J,I,K) = CF(I,J,K) 50 CONTINUE 60 CONTINUE 70 CONTINUE WRITE (ITAPE) (((CF(IORDER(I),IORDER(J),K),I=1,NCHAN),J=1,NCHAN), A K=1,LAMAX) if (iflagallcf.ne.0) then deallocate(cf) iflagallcf = 0 endif END c **** parallel **** c---------------------------------------------------------------------- c ************ subroutines for parallel code ***************** c subroutine rsct(HNP1,irowsize,icolsize) c use param use basic2 use comm_interface implicit real*8 (a-h,o-z) c-------------------------------------------------------------- c sets up and diagonalizes the hamiltonian matrix c and determines the surface amplitudes from the eigenvectors. c-------------------------------------------------------------- real*4 timei,timef,time parameter (mxn21=mznr2+1) parameter (zero=0.0d0) allocatable norder(:) common /basic1/bsto,ra,nelc,nz,lrang2,nrang2,maxnhf(mzlr2), a maxnlg(mzlr2),maxnc(mzlr2),lrang1 common /basin/eigens(mznr2,mzlr2),ends(mxn21,mzlr2),delta,eta common /cases/more,mskip,ipwinit,ipwfinal,ipolph,inast,ioneone common /cupint/mnp1,nconhp,nchan,N2HDAT,NTARSYM COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 common /inform/iread,iwrite,ipunch COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, A NBUG9 common /rcases/not1,not2,not3,nprog,ncfgp COMMON /NRBMNP/MNP2MX,IPRCENT,iprint COMMON /PART/TRACE COMMON /STATE/ENAT(MZTAR),LAT(MZTAR),ISAT(MZTAR),LPTY(MZTAR), A NAST include 'mpif.h' c common /parablock/iam,nproc dimension HNP1(irowsize,icolsize) allocatable :: value(:),vector(: , :) allocatable wmatl(:,:) real*8,allocatable :: wmatarray(:) integer :: status(MPI_STATUS_SIZE) integer, allocatable :: npointer(:,:) integer :: kcount, nksize(0:nproc-1) integer :: itagwmat,itagpoint1,itagpoint2 parameter (itagwmat=42,itagpoint1=43,itagpoint2=44) common /matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest X ,NDIV c if(iam.eq.0)write (iwrite,3000) if (not1.gt.nrang2) not1 = nrang2 noterm = not1 c.......size of the original hamiltonian matrix mnp2o=mnp1. mnp2o = noterm*nchan + ncfgp C MNP2 IS THE REDUCED NUMBER OF POLES OF ANY PARTITIONED R-MATRIX MNP2 = (NOTERM*IPRCENT)/100 MNP2 = MNP2*NCHAN + NCFGP c.......diagonalize the hamiltonian matrix and store its eigenvalues in c.......the VALUE array. The eigenvectors are distributed in the local c......."vector" arrays. c.......Then calculate the surface amplitudes. if(iam.eq.0)then IF(MNP2O.NE.MNP2)WRITE(IWRITE,3025)MNP2O,MNP2 write (iwrite,3010) c write (iwrite,3020) mnp2 endif c.......allocate eigenvalues and eigenvectors arrays c....... eigenvalues are global, eigenvectors local c....... partitioned still requires original size for value allocate(value(mnp2o),vector(irowsize,icolsize),stat=ierr) if(ierr.ne.0)stop 'allocation fails for eigenvectors' if (iam.eq.0) then if(mnp2.eq.mnp2o)then if(iprint.ge.0) X write(0,*) 'Symmetry:',mskip,' size=:',mnp2 write(35,*) 'Symmetry:',mskip,' size=:',mnp2 else if(iprint.ge.0) X write(0,*) 'Symmetry:',mskip,' size=:',mnp2o,mnp2 write(35,*) 'Symmetry:',mskip,' size=:',mnp2o,mnp2 endif endif call cpu_time(timei) if(mnp2o.gt.mnp2)then !partitioned R-matrix c Use PDSYEVX call DIAG_PDSYEVX(HNP1,irowsize,icolsize,vector,value,mnp2o,mnp2) else !full R-matrix c Use PDSYEVD call DIAG_PDSYEVD(HNP1,irowsize,icolsize,vector,value,mnp2) endif call cpu_time(timef) if (iam.eq.0) then time = timef-timei if(iprint.ge.0)write(0,2998) mskip, time/60. 2998 format(20x,'Symmetry:',i4,4x,'diagonalization time=: ' x ,1pg10.3,' min.') write(35,2999) time/60. 2999 format(2x,'diagonalization time=: ',1pg10.3,' min.') endif c C write(iwrite,*)'right after diag' c if(iam.eq.0)write(iwrite,*)(value(i),i=1,mnp2) c call flush(iwrite) c....... initialization of local surface amplitudes c idivsurface=mnp2/ndiv do 3333 III=ndiv,1,-1 iupper=mnp2-(ndiv-III)*idivsurface ilow=iupper-idivsurface iupper=ilow+idivsurface-1 if(III.eq.1)ilow=1 if(III.eq.ndiv)iupper=mnp2 C if(iam.eq.0)write(0,*)ilow,iupper allocate(wmatl(nchan,ilow:iupper),wmatarray(irowsize*icolsize), * npointer(irowsize*icolsize,2),stat=ierr) if (ierr.ne.0) stop * 'allocation fails for wmatl, wmatarray and npointer' do I=ilow,iupper do J=1,NCHAN wmatl(J,I)=zero enddo enddo do I=1,irowsize*icolsize wmatarray(I)=zero npointer(I,1)=0 npointer(I,2)=0 enddo call cpu_time(timei2) c.......determine the local surface amplitudes, wmatl(channel,eigenvalue) isizemax = noterm*nchan do icol=1,icolsize do irow=1,irowsize c............... identifying irow,icol --> I,J call demapping(irow,icol,I,J) c............... wmat is constructed only from isizemax terms C............... AND MNP2.LE.MMNP2O E-VECTORS if (I.LE.isizemax.AND.J.LE.MNP2) THEN c....... wmatl(channel,eigvn) = sum{ vector(icont,eigvn)*ends(icont,channel)} jchan = (I-1)/noterm + 1 icont = mod(I,noterm) if (icont.eq.0) icont=noterm l = l2p(jchan)+1 if((J.ge.ilow).and.(J.le.iupper))then wmatl(jchan,J) = wmatl(jchan,J) + + vector(irow,icol)*ends(icont,l) endif ENDIF enddo enddo c.......Store non-zero elements, and their co-ords, in buffers to be sent to c.......the first processor kcount=0 do J=1,mnp2 if((J.ge.ilow).and.(J.le.iupper))then do jchan=1,nchan if(wmatl(jchan,J).NE.zero)then kcount=kcount+1 wmatarray(KCOUNT) = wmatl(jchan,J) npointer(KCOUNT,1)=jchan npointer(KCOUNT,2)=J endif enddo endif enddo C7778 deallocate(vector) C write(iwrite,*)'after local surface amps' c call flush_(iwrite) c....... initialization of global surface amplitudes if (nproctest.ne.0) go to 3550 c....... Sum the local wmatl to the global wmat: c....... First, gather all values of kcount to the first processor CALL MPI_GATHER(KCOUNT,1,MPI_INTEGER,nksize(0),1,MPI_INTEGER,0, * new_comm,ierr) c....... Send all buffers to the first processor. Must make sure the c....... mpi tag on messages are different, especially the 2 integer sends. c....... Tag values are set in parameter statement above. if (iam.ne.0) then CALL MPI_SEND(wmatarray,KCOUNT,MPI_REAL8,0,itagwmat, * new_comm,ierr) if (ierr.ne.0) write(0,*)'MPI_SEND: iam, ierr = ', iam, ierr CALL MPI_SEND(npointer(1,1),KCOUNT,MPI_INTEGER,0,itagpoint1, * new_comm,ierr) if (ierr.ne.0) write(0,*)'MPI_SEND: iam, ierr = ', iam, ierr CALL MPI_SEND(npointer(1,2),KCOUNT,MPI_INTEGER,0,itagpoint2, * new_comm,ierr) if (ierr.ne.0) write(0,*)'MPI_SEND: iam, ierr = ', iam, ierr else c....... Receive the buffers from other processors, and add to global wmat do J = 1, nproc-1 CALL MPI_RECV(wmatarray,nksize(J),MPI_REAL8,J,itagwmat, * new_comm,status,ierr) if (ierr.ne.0) write(0,*)'MPI_RECV: iam, ierr = ', iam, ierr CALL MPI_RECV(npointer(1,1),nksize(J),MPI_INTEGER,J, * itagpoint1,new_comm,status,ierr) if (ierr.ne.0) write(0,*)'MPI_RECV: iam, ierr = ', iam, ierr CALL MPI_RECV(npointer(1,2),nksize(J),MPI_INTEGER,J, * itagpoint2,new_comm,status,ierr) if (ierr.ne.0) write(0,*)'MPI_RECV: iam, ierr = ', iam, ierr do L=1,nksize(J) if(((npointer(L,2)).ge.ilow).and.((npointer(L,2)).le.iupper))then wmatl(npointer(L,1),npointer(L,2)) = * wmatl(npointer(L,1),npointer(L,2)) + * wmatarray(L) endif enddo enddo endif c call mpi_barrier(new_comm,ierr) ! don't remove c c....... sum the local wmatl to the global wmat ! call mpi_reduce(wmatl,wmat,nchan*mnp2,mpi_real8,mpi_sum,0, ! + mpi_comm_world,ierr) call cpu_time(timef2) if (iam.eq.0) then time = timef2-timei2 c write(0,2998) time/60. c write(35,2998) time/60. c2998 format(2x,'global wmat sum time=:',1pg10.3,' min.') endif deallocate(wmatarray,npointer) c write(iwrite,*)'after global surface amps' c call flush_(iwrite) if (iam.ne.0) go to 3550 c.......determine EZERO for partitioned if(ndiv.eq.1)then IF(MNP2O.NE.MNP2)THEN DO I=1,MNP2 TRACE=TRACE-VALUE(I) ENDDO TRACE=TRACE/(MNP2O-MNP2) WRITE(ITAPE3)TRACE write(iwrite,3045)(TRACE-VALUE(1))*2 ENDIF endif c.......write eigenvalues to itape3 if(iupper.eq.mnp2)then write (itape3) (value(i),i=mnp2,1,-1) IF (NBUG5.NE.0 .OR. NBUG6.NE.0 .OR. NBUG8.NE.0) THEN ELOW = ENAT(1)*2.0D0 WRITE (IWRITE,3050) ELOW,(2*VALUE(I)-ELOW,I=MNP2,1,-1) ENDIF endif allocate(norder(nchan),stat=ierr) if(ierr.ne.0)stop 'allocation fails for norder' call ORDER(et,norder,nchan,-1) c c....... skip writing for dimension test mode: if (nproctest.ne.0) go to 3500 c.......write surface amplitudes to itape3 if(ndiv.eq.1)then write (itape3)((wmatl(norder(k),i),k=1,nchan),i=iupper,ilow,-1) else write (itape3)ilow,iupper,ndiv write (itape3)((wmatl(norder(k),i),k=1,nchan),i=iupper,ilow,-1) call flush(itape3) endif c write(5000,*)(k,i,(wmatl(norder(k),i),k=1,nchan),i=iupper,ilow,-1) c write (iwrite,*)wmatl(norder(1),1),wmatl(norder(nchan),mnp2) 3500 deallocate(norder) 3550 deallocate(wmatl) call mpi_barrier(new_comm,ierr) 3333 continue deallocate(value,vector) call mpi_barrier(new_comm,ierr) c write(iwrite,*)'end of RSCT' c call flush_(iwrite) c 3000 format (//35x,'SUBROUTINE RSCT'/35x,15 ('-')) 3010 format (/' DIAGONALIZING THE HAMILTONIAN MATRIX') 3020 format (/' HAMILTONIAN MATRIX - SIZE',i9/) 3025 FORMAT (/' HAMILTONIAN MATRIX SIZE REDUCED FROM',I7,' TO',I7/) 3040 format (/' EIGENVALUE/2RYD =',f14.6) 3045 FORMAT (/' PARTITIONED R-MATRIX, EZERO=',F10.3,' RYD'/) 3050 FORMAT (/' ENERGY LEVELS WITH RESPECT TO THE GROUND STATE (AT', A F12.5,' RYDBERGS)'/ (1X,10F12.5)) return end c*********************************************************************** subroutine initctxt c use param use comm_interface implicit real*8(a-h,o-z) c.......Initialize a single BLACS context common /procdat/ictxt,myrow,mycol common /matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest X ,NDIV integer i integer, allocatable :: map(:) logical :: exdis C INQUIRE(file='DISTRIBUTE.DAT',exist=exdis) if(exdis)then OPEN(444,file='DISTRIBUTE.DAT',form='formatted',status='unknown') DO I=0,NPW-1 READ(444,*)idum1,idum2,idum3,idum4 if(icolour.eq.I)istart1=idum1 ENDDO CLOSE(444) endif C allocate(map(0:nprow*npcol-1)) do i=0,nprow*npcol-1 if(.not.exdis)then map(i) = icolour * (nprow * npcol) + i ! column-major else map(i) = istart1 + i endif end do call blacs_get(0,0,ictxt) call blacs_gridmap(ictxt,map,nprow,nprow,npcol) ! column-major CALL BLACS_GRIDINFO( ictxt, NPROW, NPCOL, MYROW, MYCOL ) deallocate(map) return end c*********************************************************************** subroutine read2r c use param implicit real*8(a-h,o-z) common /crashblock/istop c.......read "sizeH.dat" and determine the global size of H c c....... or sizeBP.dat c parameter (mxmat=mzchf*mznr2+mznc2) common /inform/iread,iwrite,ipunch common /globsize/idimH(mzslp) character a1*8,a2*9,a3*7 logical ex ex=.true. inquire (file='sizeBP.dat',exist=ex) if (ex) then c........ Breit-Pauli mode open(unit=30,file='sizeBP.dat',status='old') i=1 80 read(30,*,err=80,end=1000) nchan,nconhp,mnp1 c write(0,*)'reading nchan=',nchan,' mnp1',mnp1 if (mnp1.gt.mxmat) call recov2('read2r','mxmat ',mxmat,mnp1) if (istop.ne.0) return idimH(i) = mnp1 i=i+1 if (i.gt.mzslp) call recov2('read2r','mzslp ',mzslp,i) if (istop.ne.0) return go to 80 else c......... ICFT mode inquire (file='sizeH.dat',exist=ex) if (ex) then open(unit=30,file='sizeH.dat',status='old') i=1 100 read(30,3010,err=100,end=1000) a1,nchan,a2,nconhp,a3,mnp1 if (a1(1:8).ne.' NCHAN ='.or.a2(1:9).ne.' NCONHP ='.or. + a3(1:7).ne.' MNP1 =') go to 100 if (mnp1.gt.mxmat) call recov2('read2r','mxmat ',mxmat,mnp1) if (istop.ne.0) return idimH(i) = mnp1 cold if (nchan.eq.0) go to 100 ! as now KOUNTed if (nchan.eq.0) idimH(i) = 0 ! mnp1>0 if corr i=i+1 if (i.gt.mzslp) call recov2('read2r','mzslp ',mzslp,i) if (istop.ne.0) return go to 100 3010 FORMAT (a8,I5,a9,I7,a7,I7) endif !! ICFT endif if (.not.ex) then write(iwrite,7777) istop=80 return 7777 format (/' Neither routjk or sizeH.dat exist....stopping.') endif 1000 close(30) idimH(i)=-1 !flag eof return end c*********************************************************************** subroutine Hsize(irowsize,icolsize) c use param implicit real*8(a-h,o-z) c.......calculate the size of the local Hamiltonian include 'mpif.h' common /matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest X ,NDIV common /procdat/ictxt,myrow,mycol common /globsize/idimH(mzslp) common /cases/more,mskip,ipwinit,ipwfinal,ipolph,inast,ioneone cv3 call MPI_BCAST(NB,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) cv3 call MPI_BCAST(MB,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) cv3 call MPI_BCAST(nprow,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) cv3 call MPI_BCAST(npcol,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) cv3 call MPI_BCAST(irsrc,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) cv3 call MPI_BCAST(icsrc,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) cv3 call MPI_BCAST(itype,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) NN=idimH(mskip) c.......dimension of the local matrix HNP1 irowsize = numroc(NN,nb,myrow,irsrc,nprow) icolsize = numroc(NN,mb,mycol,icsrc,npcol) C write(0,*)'inside numroc',NN,nb,myrow,irsrc,icsrc,nprow,npcol return end c*********************************************************************** subroutine diag_pdsyevx(HNP1,irowsize,icolsize,Z,W,NN,IU) c use param use comm_interface implicit real*8(a-h,o-z) c.......Matrix diagonalization for parallel machine c.......driver for the SCAlapack subroutine pdsyevx include 'mpif.h' C common /parablock/iam,nproc common /procdat/ictxt,myrow,mycol dimension HNP1(irowsize,icolsize) dimension Z(irowsize,icolsize) ! local eigenvectors dimension W(NN) ! global eigenvalues parameter (zero=0.0d0) c.......parameters and variables used in pdsyevx parameter (rminusone=-1.0d0) parameter (abstol = rminusone) ! error tolerance of eigenvalues parameter (orfac = 1.d-6) ! reorthogonalization parameter (iclustersize=11) integer desca(9),descz(9) allocatable :: ifail(:),iwork(:),iclustr(:) allocatable :: work(:),gap(:) common /matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest X ,NDIV c.......other parameters needed in pdsyevx lda = max(irowsize,1) allocate(ifail(NN),stat=ierr) if(ierr.ne.0)stop 'allocation fails for ifail' c.......space required is such that all eigenvectors could be computed, c.......even when partitioned requires less. neig = NN NNN = max(NN,nb,2) np0 = numroc(NNN,nb,mod(irsrc,nprow),mod(irsrc,nprow),nprow) mmm=mod(icsrc,npcol) mq0 = numroc(max(neig,nb,2),nb,mmm,mmm,npcol) lwork = 5*NN + max(5*NNN,np0*mq0+2*nb*nb) + + iceil(neig,nproc)*NNN + (iclustersize-1)*NN allocate(work(lwork),stat=ierr) if(ierr.ne.0)stop 'allocation fails for work' c.......working array liwork = max(3*NN+nproc+1,4*NN,nproc,14)+2*NN allocate(iwork(liwork),stat=ierr) if(ierr.ne.0)stop 'allocation fails for iwork' c.......gap and iclustr allocate(iclustr(nproc*2),stat=ierr) if (ierr.ne.0)stop 'allocation fails for iclustr' allocate(gap(nproc),stat=ierr) if (ierr.ne.0)stop 'allocation fails for gap' if(nproctest.ne.0)go to 6666 c.......basic array descriptors call descinit(desca,NN,NN,nb,mb,irsrc,icsrc,ictxt,lda,info) call descinit(descz,NN,NN,nb,mb,irsrc,icsrc,ictxt,lda,info) call mpi_barrier(new_comm,ierr) c......compute the partial eigendecomposition in PDSYEVX c pdsyevx for SUN, SGI and IBM --- pssyevx for CRAY call pdsyevx('V','I','U',NN,HNP1,1,1,desca,zero,zero,1, + IU,abstol,m,nz,W,orfac,Z,1,1,descz,work, + lwork,iwork,liwork,ifail,iclustr,gap,info) if(m.ne.IU)then if (iam.eq.0) x write(0,*) ' pdsyevx error - not all e-values found:', m,iu write(0,*)' pdsyevx error - not all e-values found:',m,iu stop 'pdsyevx error' endif if(m.ne.nz)then if (iam.eq.0) x write(0,*) ' pdsyevx error - not all e-vectors found:', m,nz write(0,*)' pdsyevx error - not all e-vectors found:',m,nz stop 'pdsyevx error' endif if(info.ne.0) then if (iam.eq.0) write(0,*) ' pdsyevx error - info=:',info write(0,*) ' pdsyevx error - info=:',info if(mod(info/2,2).ne.0)then write(0,*)'Try increasing iclustersize, cluster info:' write(0,*)(iclustr(i),i=1,2*nproc) endif stop 'pdsyevx error' endif 6666 deallocate(ifail) deallocate(work) deallocate(iwork) deallocate(gap) deallocate(iclustr) return end c*********************************************************************** subroutine diag_pdsyevd(HNP1,irowsize,icolsize,Z,W,NN) c use param use comm_interface implicit real*8(a-h,o-z) c.......Matrix diagonalization for parallel machine c.......driver for the SCAlapack subroutine pdsyevd include 'mpif.h' C common /parablock/iam,nproc common /procdat/ictxt,myrow,mycol dimension HNP1(irowsize,icolsize) dimension Z(irowsize,icolsize) ! local eigenvectors dimension W(NN) ! global eigenvalues c.......parameters and variables used in pdsyevd integer desca(9),descz(9) integer iwork allocatable :: iwork(:) allocatable :: work(:) integer ione,itwo,ithree,isix,iseven,ieight data ione,itwo,ithree,isix,iseven,ieight/1,2,3,6,7,8/ common /matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest X ,NDIV c.......other parameters needed in pdsyevd lda = max(irowsize,ione) c.......space required to guarantee that all eigenvectors are computed np0 = numroc(NN,nb,myrow,irsrc,nprow) mq0 = numroc(NN,nb,mycol,icsrc,npcol) itril = ithree*NN + max(nb*(np0+ione),ithree*nb) lwork = itwo*NN + max(ione+isix*NN+itwo*np0*mq0,itril) allocate(work(lwork),stat=ierr) if(ierr.ne.0)stop 'allocation fails for work' c.......working array liwork = iseven*NN + ieight*npcol + itwo allocate(iwork(liwork),stat=ierr) if(ierr.ne.0)stop 'allocation fails for iwork' if(nproctest.ne.0)go to 6666 c.......basic array descriptors c write(0,*)'iam=',iam,desca,NN,NN,nb,mb,irsrc,icsrc,ictxt call descinit(desca,NN,NN,nb,mb,irsrc,icsrc,ictxt,lda,info) call descinit(descz,NN,NN,nb,mb,irsrc,icsrc,ictxt,lda,info) call mpi_barrier(new_comm,ierr) c......compute the entire eigendecomposition in PDSYEVD c pdsyevd SUN, SGI and IBM, pssyevd CRAY call pdsyevd('V','U',NN,HNP1,1,1,desca,W,Z,1,1,descz,work, + lwork,iwork,liwork,info) if(info.ne.0) then if (iam.eq.0) write(0,*) ' pdsyevd error - info=:',info stop 'pdsyevd error' endif 6666 deallocate(work) deallocate(iwork) return end c*********************************************************************** subroutine demapping(irow,icol,I,J) c use param implicit real*8(a-h,o-z) c.......demapping the local (irow,icol) in the global (I,J) position common /matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest X ,NDIV common /procdat/ictxt,myrow,mycol ix = mod(irow,nb) if (ix.eq.0) ix=nb iy = mod(icol,mb) if (iy.eq.0) iy=mb c....... l and m: local block coordinate l = (irow-ix)/nb m = (icol-iy)/mb I = (l*nprow + myrow) * nb + ix J = (m*npcol + mycol) * mb + iy return end c*********************************************************************** subroutine locH(ijgbindex,irow,icol,ipr,ipc) c use param implicit real*8(a-h,o-z) integer*8 ijgbindex c.......Mapping the global matrix into the local matrix c....... IJGBINDEX ----> (I,J) ------> (irow,icol) common /cases/more,mskip,ipwinit,ipwfinal,ipolph,inast,ioneone common /globsize/idimH(mzslp) common /matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest X ,NDIV common /procdat/ictxt,myrow,mycol c........identifying I,J from ijgbindex NN = idimH(mskip) call ntoij(NN,ijgbindex,i,j) c........mapping the GLOBAL matrix A(I,J) into LOCAL A(irow,icol) c........ipr and ipc: row and column processes call mapping(I,J,irow,icol,ipr,ipc) return end c*********************************************************************** subroutine ntoij(NN,IJ,i,j) c use param implicit none integer :: NN,I,J real*8 :: b,bb,q,two integer*8 ibig,jbig,IJ,NNbig,NNbig2 parameter(two=2.0d0) c........identifying i,j from index IJ NNbig=NN NNbig2=2*NNbig b = (NNbig2+3)/two bb=b*b q=two*(NNbig+IJ) ibig = b - sqrt(bb-q) jbig = NNbig + IJ - ibig*(NNbig2 - ibig + 1)/two i=ibig j=jbig return end c*********************************************************************** subroutine ijton(NN,IJ,i,j) c use param implicit real*8(a-h,o-z) integer*8 IJ,NNbig,ibig,jbig c.......identifying index IJ from i,j ibig=i jbig=j NNbig=NN c IJ = (NNbig*(NNbig+1))/2-((NNbig-ibig+1)*(NNbig-ibig+2))/2+jbig-ibig+1 IJ = JBIG+(IBIG-1)*(2*NNbig-IBIG)/2 return end c*********************************************************************** subroutine mapping(I,J,irow,icol,ipr,ipc) c use param implicit real*8(a-h,o-z) c.......mapping the GLOBAL matrix A(I,J) into LOCAL A(irow,icol) common /matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest X ,NDIV common /procdat/ictxt,myrow,mycol c IM=I-1 JM=J-1 c..........ipr and ipc: row and column processes ipr = mod(irsrc + IM/MB , nprow) ipc = mod(icsrc + JM/NB , npcol) c.......column identification (J --> icol) mblock = JM/(npcol*NB) iypos = mod(JM,NB) + 1 icol = mblock*NB + iypos c.......row identification (I --> irow) lblock = IM/(nprow*MB) ixpos = mod(IM,MB) + 1 irow = lblock*MB + ixpos return end c*********************************************************************** c **** parallel **** subroutine nochan() use param use comm_interface implicit real*8(a-h,o-z) C C Removes empty STG2HXXX.DAT files for the cases C where no channels couple for a particular symmetry C C C character a1*8,a2*9,a3*7 integer :: imiss(0:MZSLP)=0,nchana(0:MZSLP)=0, X nconhpa(0:MZSLP)=0,mnp1b(0:MZSLP)=0 CHARACTER :: NUM(0:9)*1 CHARACTER*12 :: stg2a,stg2b DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ C C C open(631,file='sizeH.dat',form='formatted',status='old') if(masteriam.eq.0)then open(630,file='sizeH.datb',form='formatted',status='unknown') endif C i=0 100 read(631,3010,err=100,end=1000)a1,nchana(i),a2,nconhpa(i),a3, X mnp1b(i) if(masteriam.eq.0.and.nchana(i).ne.0)then write(630,3010)a1,nchana(i),a2,nconhpa(i),a3,mnp1b(i) endif if (nchana(i).eq.0) then if(masteriam.eq.0)write(0,*)'Symmetry ',i,': has 0',' channels' imiss(i)=1 endif i=i+1 goto 100 1000 continue C KOUNT=I C IFLAG=0 C IFLAG=MAXVAL(IMISS) C IF(IFLAG.eq.0)THEN CLOSE(630) CLOSE(631) RETURN ENDIF C IF(MASTERIAM.EQ.0)WRITE(0,*) X 'PARTIAL WAVES WITH NO COUPLED CHANNELS DISCOVERED' C IF(npw_per_subworld.eq.1)then !TBD should be IONEONE=1 now... if(masteriam.eq.0)then write(0,*)'Will correct sizeH.dat and STG2HXXX.DAT ' write(0,*)'Rerun with number of processors accordingly reduced' endif else if(masteriam.eq.0)then write(0,*)'setting npw_per_subworld= wc sizeH.dat' write(0,*)'is the safest approach but will continue ' endif close(631) return ENDIF CLOSE(631) C IF (MASTERIAM.EQ.0)THEN C 1200 CONTINUE DO J=0,kount-1 if(imiss(J).eq.1)then do N=J+1,kount i1=(N-1)/100 i2=((N-1)-100*((N-1)/100))/10 i3=(N-1)-(100*((N-1)/100))-i2*10 stg2a='STG2H'//num(i1)//num(i2)//num(i3)//'.DAT' C i4=N/100 i5=(N-100*(N/100))/10 i6=N-(100*(N/100))-i5*10 stg2b='STG2H'//num(i4)//num(i5)//num(i6)//'.DAT' C open(632,file='script',form='formatted',status='unknown') rewind(632) write(632,*)"mv ",stg2b," ",stg2a close(632) C call system("chmod u+x script") call system("./script") C imiss(N-1)=imiss(N) C enddo endif kount=kount-1 ENDDO C close(630) C open(632,file='script',form='formatted',status='unknown') rewind(632) write(632,*)"mv sizeH.datb sizeH.dat" close(632) C call system("chmod u+x script") call system("./script") C 3010 FORMAT (a8,I5,a9,I7,a7,I7) ENDIF C call comm_barrier() call comm_finalize() stop end c*********************************************************************** subroutine ntoijandback(NN,IJ,iiij) c use param implicit real*8(a-h,o-z) integer*8 ibig,jbig,IJ,NNbig,iiij parameter(two=2.0d0) c........identifying i,j from index IJ NNbig=NN b = (2*NNbig+3)/two bb=b*b q=two*(NNbig+IJ) ibig = b - sqrt(bb-q) jbig = NNbig + IJ - ibig*(2*NNbig - ibig + 1)/two c i=ibig c j=jbig IIIJ = JBIG+(IBIG-1)*(2*NNbig-IBIG)/2 return end c***********************************************************************