DDotMPI and DGEMVMPI implementation

.gitignore

@@ -1,4 +1,5 @@
 *.o
 *.mod
+*.swp
 lib/*
 bin/*

Changes_cychen.txt

@@ -0,0 +1,17 @@
+00. src/tool: Makefile, add "chmod a+x $(GRASP)/bin/rsave"
+01. jj2lsj90: In jj2lsj_code.f90, lines 220 and 222 are changed to deal with both parities present. 
+#02. jj2lsj90: In jj2lsj_code.f90, lines around 850, and lines around 2435,
+#modifications for transformation of configurations with empty K-shell, such as
+#2s2p, 2lnl, as done in A&A 592, A141(2016). 
+03. rcsfzerofirst90: In RCSFzerofirst.f90, the number of CSFs in zero-order space
+are output. The informations could be used conveniently in RCSF and RCI calculations.
+04. rmcdhf90: In newco.f90, weighted average energy is output to monitor the energy convergence.
+05. rangular90_mpi: getinf.90 modified for iccut=1
+06. rci90_mpi: setham_gg.f90, only myid=0 reports how far the calculation has proceeded. 
+07. libdvd90, dvdson.f90: make lines 1119 and 1120 work for convergence
+monitor. The output messages are very useful for extremely large-scaled calculations. 
+08. libdvd90, dvdson.f90: function TSTSEL is modified, so that dvdson performs at least TWO
+iteration calculations. Or, error results could be obtained in some cases (in
+O-like isoelectronic sequence?). 
+09. dvdson--mpi
+10. Modifcations for searching the targeted eigenpairs one by one

make_environment_gfortran

@@ -9,13 +9,20 @@
 #
 # Installation requirements:
 # - Lapack, Blas and MPI libraries have to be installed and properly linked - e.g. add them to LD_LIBRARY_PATH.
+mkdir lib bin > /dev/null 2>&1 
+echo "cp ../BLAS-3.8.0/blas_LINUX.a lib/libblas.a"
+cp ../BLAS-3.8.0/blas_LINUX.a lib/libblas.a
+echo "cp ../lapack-3.8.0/liblapack.a lib/"
+cp ../lapack-3.8.0/liblapack.a lib/
+
 # - The Fortran compiler of choice and the MPI wrapper (as specified by FC and FC_MPI below) have to be on your PATH.
 #
 # -------------------------------------------------------------------------------------------------------------------
 # Set up main flags
 # -------------------------------------------------------------------------------------------------------------------
 export   FC=gfortran                                           # Fortran compiler
-export   FC_FLAGS="-O2 -fno-automatic "                        # Serial code compiler flags
+#export   FC_FLAGS="-O2 -fno-automatic"                        # Serial code compiler flags
+export   FC_FLAGS="-O2 -fno-automatic -fconvert=big-endian"    # If little-endian used, comment out this line, use the above one
 export   FC_LD=" "                                             # Serial linker flags
 export   GRASP="${PWD}"                                        # Location of the 2018 root directory
 export   LAPACK_LIBS="-llapack -lblas"                         # Lapack libraries
@@ -26,4 +33,4 @@ export   FC_MPI="mpifort"                                      # MPI
 export   FC_MPIFLAGS="${FC_FLAGS}"                             # Parallel code compiler flags
 export   FC_MPILD=${FC_LD}                                     # Serial linker flags
 # -------------------------------------------------------------------------------------------------------------------
-export   MPI_TMP="${HOME}/grasp_mpi_tmp"                       # Location for temporary files
+export   MPI_TMP="${HOME}/tmp"                       # Location for temporary files

src/appl/Makefile

@@ -1,7 +1,8 @@
-SUBDIR =  HF jj2lsj90 jjgen90 rangular90  rbiotransform90  rci90  rcsfgenerate90  \
-          rcsfinteract90 rcsfzerofirst90 rhfs90 rmcdhf90 rnucleus90 rtransition90  rwfnestimate90 sms90  \
+SUBDIR =  HF jj2lsj90 jjgen90 rangular90  rbiotransform90  rcsfgenerate90  \
+          rcsfinteract90 rcsfzerofirst90 rhfs90 rnucleus90 rtransition90  rwfnestimate90 sms90  \
 	  rangular90_mpi rbiotransform90_mpi rci90_mpi rmcdhf90_mpi rtransition90_mpi
-
+          # rci90 rmcdhf90 
+
 TARGETS	 = install
 
 $(TARGETS):

src/appl/rangular90_mpi/getinf.f90

@@ -60,8 +60,9 @@ SUBROUTINE GETINF
             do i = 1,nblock
               write(istde,*) 'Give ICCUT for block',i
     1         READ *, ICCUT(i)
-              IF ((ICCUT(i) <= 1).OR.(ICCUT(i) >= ncfblk(i))) THEN
-                WRITE (istde,*) 'GETINF: ICCUT must be greater than 1',&
+!cychen              IF ((ICCUT(i) <= 1).OR.(ICCUT(i) >= ncfblk(i))) THEN
+              IF ((ICCUT(i) < 1).OR.(ICCUT(i) >= ncfblk(i))) THEN
+                WRITE (istde,*) 'GETINF: ICCUT must be greater than 0',&
                                 ' and less than ',ncfblk(i)
                 WRITE (istde,*) ' please reenter ICCUT:'
                 GOTO 1

src/appl/rci90/maneig.f90

@@ -107,6 +107,10 @@ SUBROUTINE MANEIG(IATJPO, IASPAR)
       REAL(DOUBLE), DIMENSION(:), pointer :: w, z, work, diag
       INTEGER, DIMENSION(:), pointer :: iwork, ifail, jwork
 !-----------------------------------------------------------------------
+!CYC: Search the targeted eigenpairs one by one
+      INTEGER IRESTART_GDVD
+      IRESTART_GDVD = 0
+
       ABSTOL = 2*DLAMCH('S')
       MYID = 0
       NPROCS = 1
@@ -436,7 +440,7 @@ SUBROUTINE MANEIG(IATJPO, IASPAR)
 
                !NIV = 0   ! Why equal 0 ???
                !WRITE (6, *) ' Calling gdvd(spodmv,...'
-               CALL GDVD (SPODMV, NCF, LIM, DIAG, ILOW, IHIGH, JWORK, NIV, &
+               CALL GDVD (SPODMV, IRESTART_GDVD, NCF, LIM, DIAG, ILOW, IHIGH, JWORK, NIV, &
                   MBLOCK, CRITE, CRITC, CRITR, ORTHO, MAXITR, WORK, LWORK, &
                   IWORK, LIWORK, HIEND, NLOOPS, NMV, IERR)
 
@@ -448,7 +452,7 @@ SUBROUTINE MANEIG(IATJPO, IASPAR)
                !WRITE (*, *) NCF, NIV, (WORK(I),I=NCF*NIV + 1,NCF*NIV + NIV)
                !WRITE (*, *) LIM, ILOW, IHIGH, MBLOCK, MAXITR, LWORK, LIWORK
                !WRITE (*, *) IERR
-               CALL GDVD (SPICMV2, NCF, LIM, DIAG, ILOW, IHIGH, JWORK, NIV, &
+               CALL GDVD (SPICMV2, IRESTART_GDVD, NCF, LIM, DIAG, ILOW, IHIGH, JWORK, NIV, &
                   MBLOCK, CRITE, CRITC, CRITR, ORTHO, MAXITR, WORK, LWORK, &
                   IWORK, LIWORK, HIEND, NLOOPS, NMV, IERR)
                !WRITE (*, *) 'after gdvd...'
@@ -465,7 +469,7 @@ SUBROUTINE MANEIG(IATJPO, IASPAR)
                WRITE (6, *) ' Dense - Memory,  iniestdm'
 !            CALL INIESTDM (1000,NCF,NIV,WORK,EMT)
                CALL INIESTDM (2000, NCF, NIV, WORK, EMT)
-               CALL GDVD (DNICMV, NCF, LIM, DIAG, ILOW, IHIGH, JWORK, NIV, &
+               CALL GDVD (DNICMV, IRESTART_GDVD, NCF, LIM, DIAG, ILOW, IHIGH, JWORK, NIV, &
                   MBLOCK, CRITE, CRITC, CRITR, ORTHO, MAXITR, WORK, LWORK, &
                   IWORK, LIWORK, HIEND, NLOOPS, NMV, IERR)
                CALL DALLOC (EMT, 'EMT', 'MANEIG')

src/appl/rci90_mpi/maneigmpi.f90

@@ -132,8 +132,12 @@ SUBROUTINE MANEIG(IATJPO, IASPAR, NELMNT_a)
       CHARACTER(LEN=8) :: CNUM
       REAL(DOUBLE), DIMENSION(:), pointer :: w, z, work, diag
       INTEGER, DIMENSION(:), pointer :: iwork, ifail, jwork
-!-----------------------------------------------
-!
+
+!-----------------------------------------------------------------------
+! CYC Modifications for step by step diagonalization
+      INTEGER IRESTART_GDVD, ISTEP, NSTEP, IITMP, NMV0
+      INTEGER TIME0, TIME1, TIME2
+! CYC
 !-----------------------------------------------------------------------
       ABSTOL = 2*DLAMCH('S')
       IF (MYID == 0) WRITE (6, *) 'Calling maneig...'
@@ -494,7 +498,7 @@ SUBROUTINE MANEIG(IATJPO, IASPAR, NELMNT_a)
               print *, 'Returned from iniestsd '
               if (ncf.gt. IOLPCK) then
                 print *, 'Calling GDVD'
-                CALL GDVD (SPODMV,NCF,LIM,DIAG,ILOW,IHIGH,            &
+                CALL GDVD (SPODMV,IRESTART_GDVD,NCF,LIM,DIAG,ILOW,IHIGH,            &
                   JWORK,NIV,MBLOCK,CRITE,CRITC, CRITR,ORTHO,MAXITR,   &
                   WORK,LWORK,IWORK,LIWORK,HIEND,NLOOPS,               &
                   NMV,IERR)
@@ -505,10 +509,70 @@ SUBROUTINE MANEIG(IATJPO, IASPAR, NELMNT_a)
                IF (myid .EQ. 0) print *, ' Sparse - Memory, iniestmpi'
                CALL iniestmpi (IOLPCK, NCF,NIV,WORK,EMT,IENDC,IROW)
                if(ncf.gt. IOLPCK) then
-                 CALL GDVD (SPICMVmpi,NCF,LIM,DIAG,ILOW,IHIGH,        &
+!-----------------------------------------------------------------------
+! CYC: Try to diagonalize step by step
+                IRESTART_GDVD = 0
+                TIME0 = MPI_WTIME()
+! More tests should be performed for istep, it could be set as 1
+                ISTEP = 5 
+                NMV = 0 
+                IF (NCF.GT.2.0E5 .AND. NVEX.GT.ISTEP) THEN
+                 NSTEP = NVEX / ISTEP
+                 IF (MOD(NVEX,ISTEP).LE.ISTEP/2) NSTEP = NSTEP - 1
+                 DO IITMP = 1, NSTEP
+! NIV could be changed to nvecmx+1 in dvdson:initdvd procedure
+                  NIV = NVEX
+                  ILOW = 1
+                  IHIGH = IITMP * ISTEP
+                  IF (IITMP .EQ. 1 ) THEN
+                    IRESTART_GDVD = 0
+                  ELSE
+                    IRESTART_GDVD = 1
+                  ENDIF
+                  NMV0 = NMV
+                  TIME1 = MPI_WTIME()
+
+                  CALL GDVD (SPICMVmpi,IRESTART_GDVD,NCF,LIM,DIAG,ILOW,IHIGH, &
                   JWORK,NIV,MBLOCK,CRITE,CRITC, CRITR,ORTHO,MAXITR,   &
                   WORK,LWORK,IWORK,LIWORK,HIEND,NLOOPS,               &
                   NMV,IERR)
+
+                  TIME2 = MPI_WTIME()
+                  IF (MYID .EQ. 0) THEN
+                    WRITE (*, *)' IITMP= ', IITMP,  &
+                                ' Time (s)=', TIME2 - TIME1
+                    WRITE (*,*) ' ', NMV - NMV0, &
+                      'Matrix-vector multiplies for this loop.'  
+                    WRITE (*,*) ' ', NMV,  &
+                      ' Total matrix-vector multiplies by now.'
+                  ENDIF
+                 ENDDO
+                ENDIF
+
+! Search the remaining (NVEX - ISTEP * NSTEP) EIGENPAIRS
+                NIV = NVEX
+                ILOW = 1
+                IHIGH = NVEX
+                NMV0 = NMV  
+                TIME1 = MPI_WTIME()
+
+                CALL GDVD(SPICMVmpi,IRESTART_GDVD,NCF,LIM,DIAG,ILOW,IHIGH, &
+                JWORK,NIV,MBLOCK,CRITE,CRITC, CRITR,ORTHO,MAXITR,   &
+                WORK,LWORK,IWORK,LIWORK,HIEND,NLOOPS,               &
+                NMV,IERR)
+
+                TIME2 = MPI_WTIME()
+                IF (MYID .EQ. 0) THEN
+                  WRITE (*, *)' IITMP= ', IITMP,  &
+                              ' Time (s)=', TIME2 - TIME1 
+                  WRITE (*,*) ' ', NMV - NMV0, &
+                    'Matrix-vector multiplies for this loop.'  
+                  WRITE (*,*) ' ', NMV,  &
+                    ' Total matrix-vector multiplies by now.'
+                  WRITE (*, *)' Total time (m)=', (TIME2 - TIME0) / 60.0
+                ENDIF
+! CYC
+!-----------------------------------------------------------------------
                end if
 
                CALL DALLOC (EMT, 'EMT', 'MANEIG')
@@ -520,7 +584,7 @@ SUBROUTINE MANEIG(IATJPO, IASPAR, NELMNT_a)
                IF (myid .EQ. 0) print *, ' Dense - Memory, iniestdm'
               CALL INIESTDM (IOLPCK,NCF,NIV,WORK,EMT)
               if (ncf.gt. IOLPCK) then
-                CALL GDVD (DNICMV,NCF,LIM,DIAG,ILOW,IHIGH,            &
+                CALL GDVD (DNICMV,IRESTART_GDVD,NCF,LIM,DIAG,ILOW,IHIGH,            &
                     JWORK,NIV,MBLOCK,CRITE,CRITC, CRITR,ORTHO,MAXITR, &
                     WORK,LWORK,IWORK,LIWORK,HIEND,NLOOPS,             &
                     NMV,IERR)

src/appl/rci90_mpi/setham_gg.f90

@@ -398,6 +398,8 @@ SUBROUTINE SETHAM (myid, nprocs, jblock, ELSTO,ICSTRT, nelmntt, &
 !cjb MPI progress begin --------------------------------------------------
 !
 !cjb just started
+!cychen: report only from myid.eq.0
+      if (myid .eq. 0 ) then
        IF (IC .LE. nprocs) then
             PRINT '(A6,I10,A9,I10,A9,I5,A8,I5)', 'Start ', ic,       &
                     ' nnonz = ', nelc,                             &
@@ -410,19 +412,19 @@ SUBROUTINE SETHAM (myid, nprocs, jblock, ELSTO,ICSTRT, nelmntt, &
        IF (IC .GT. nprocs .and. IC .LT. NCF-nprocs .and. &
            MOD (IC-1,100*nprocs) .EQ. myid) then
          if (myid .eq. 0) then
-            PRINT '(A5,I10,A9,I10,A9,I5,A8,I5)', 'Done ', ic,       &
+            PRINT '(A6,I10,A9,I10,A9,I5,A8,I5)', 'Done ', ic,       &
                     ' nnonz = ', nelc,                             &
                     ' block = ', jblock,' myid = ', myid
         flush output_unit
         flush error_unit
          else
-            PRINT '(A5,I10,A9,I10,A9,I5,A8,I5)', 'Row  ', ic,       &
+            PRINT '(A6,I10,A9,I10,A9,I5,A8,I5)', 'Row  ', ic,       &
                     ' nnonz = ', nelc,                             &
                     ' block = ', jblock,' myid = ', myid
         flush output_unit
         flush error_unit
          endif
-       endif
+      endif
 !
 !cjb almost finished
        IF (IC .GT. NCF-nprocs) then
@@ -433,6 +435,7 @@ SUBROUTINE SETHAM (myid, nprocs, jblock, ELSTO,ICSTRT, nelmntt, &
         flush output_unit
         flush error_unit
        endif
+      endif
 !cjb MPI progress end ----------------------------------------------------
 !
 !

src/appl/rcsfzerofirst90/RCSFzerofirst.f90

@@ -44,19 +44,23 @@ PROGRAM RCSFzerofirst
       print *, ""
       NBLOCK = 0
       CALL SET_CSF_ZFlist
+!cychen: output the zero-order space for re-use in mcp, rci
+      open(301,file='icut',form='formatted',status='unknown')
       WRITE (6, *) "  Block    Zero-order Space   Complete Space"
       DO
          CALL LODCSL_Zero (NEXT_BLOCK)
          CALL LODCSL_Part (CSF_Number)
          WRITE (6,'(3X,I2,6X,I14,3X,I17)')                            &
                                      NBLOCK,NCFBLK(NBLOCK),CSF_Number-1
+         write(301,*)NCFBLK(NBLOCK)
          deallocate (Found)
          deallocate (C_shell)
          deallocate (C_quant)
          deallocate (C_coupl)
          IF(.NOT. NEXT_BLOCK) EXIT
          WRITE(22,'(A2)') ' *'
       END DO
+      close(301)
       call stoptime (ncount1, 'RCSFzerofirst')
       STOP
       END PROGRAM RCSFzerofirst

src/appl/rmcdhf90/maneig.f90

@@ -78,6 +78,10 @@ SUBROUTINE MANEIG(dvdfirst, LPRINT, JBLOCK, &
       REAL(DOUBLE), DIMENSION(:), POINTER :: WORK, DIAG
       LOGICAL :: HIEND
 !-----------------------------------------------
+!CYC: Search the targeted eigenpairs one by one
+      INTEGER IRESTART_GDVD
+      IRESTART_GDVD = 0
+!
 
       !PRINT *, 'maneig ...'
 
@@ -159,7 +163,8 @@ SUBROUTINE MANEIG(dvdfirst, LPRINT, JBLOCK, &
 
          CALL ALLOC (IWORK, LIWORK, 'IWORK', 'MANEIG')
          CALL ALLOC (JWORK, LIM, 'JWORK', 'MANEIG')
-         CALL GDVD (SPICMV2, NCF, LIM, DIAG, ILOW, IHIGH, JWORK, NIV, MBLOCK, &
+         !CALL GDVD (SPICMV2, NCF, LIM, DIAG, ILOW, IHIGH, JWORK, NIV, MBLOCK, &
+         CALL GDVD (SPICMV2, IRESTART_GDVD, NCF, LIM, DIAG, ILOW, IHIGH, JWORK, NIV, MBLOCK, &
             CRITE, CRITC, CRITR, ORTHO, MAXITR, WORK, LWORK, IWORK, LIWORK, &
             HIEND, NLOOPS, NMV, IERR)
 

src/appl/rmcdhf90/newco.f90

@@ -114,6 +114,8 @@ SUBROUTINE NEWCO(SUM)
 !   Write out average energy
 !
       IF (NCMIN > 1) WRITE (*, 304) SUM
+!cychen, output sum to the terminal for convergence monitor.
+      WRITE (0,'(A25,1PD18.10)') 'Weighted average energy: ', SUM
 !
 !   Write out generalized occupation numbers
 !

src/appl/rmcdhf90_mpi/maneigmpi.f90

@@ -78,8 +78,12 @@ SUBROUTINE MANEIGmpi(dvdfirst, LPRINT, JBLOCK, &
       REAL(DOUBLE) :: PNWORK, CRITE, CRITC, CRITR, ORTHO, AMAX, WA, DNFAC
       REAL(DOUBLE), DIMENSION(:), POINTER :: WORK, DIAG, atmp
       LOGICAL :: HIEND
-!-----------------------------------------------
 
+!-----------------------------------------------
+!CYC: Search the targeted eigenpairs one by one
+      INTEGER IRESTART_GDVD
+      IRESTART_GDVD = 0
+!-----------------------------------------------
       !PRINT *, 'maneig ...'
 
 !      ...spicmvmpi needs this COMMON /WCHBLK/JBLOCKK
@@ -174,7 +178,8 @@ SUBROUTINE MANEIGmpi(dvdfirst, LPRINT, JBLOCK, &
          CALL ALLOC (IWORK, LIWORK, 'IWORK', 'MANEIGmpi')
          CALL ALLOC (JWORK, LIM, 'JWORK', 'MANEIGmpi')
          if (ncf.gt.1000) then
-            CALL GDVD (SPICMVMPI,NCF,LIM,DIAG,ILOW,IHIGH,JWORK,NIV,MBLOCK, &
+            !CALL GDVD (SPICMVMPI,NCF,LIM,DIAG,ILOW,IHIGH,JWORK,NIV,MBLOCK, &
+            CALL GDVD (SPICMVMPI,IRESTART_GDVD,NCF,LIM,DIAG,ILOW,IHIGH,JWORK,NIV,MBLOCK, &
             CRITE, CRITC, CRITR, ORTHO, MAXITR, WORK, LWORK, IWORK, LIWORK, &
             HIEND, NLOOPS, NMV, IERR)
          end if