diff --git a/arch/CRAY-T3E.pdbg b/arch/CRAY-T3E.pdbg index 41eb2fb613..67d591fe44 100644 --- a/arch/CRAY-T3E.pdbg +++ b/arch/CRAY-T3E.pdbg @@ -5,26 +5,21 @@ FC = f90 -f free FC_fixed = f90 -f fixed LD = f90 AR = ar -r -CPPFLAGS = -C -D__T3E -D__FFTSG -D__parallel -P\ +DFLAGS = -D__T3E -D__FFTSG -D__parallel\ -Ddcopy=scopy -Ddgbsv=sgbsv -Ddgecon=sgecon -Ddgemm=sgemm\ -Ddgemv=sgemv -Ddger=sger -Ddgerfs=sgerfs -Ddgetrf=sgetrf\ -Ddgetri=sgetri -Ddgetrs=sgetrs -Ddlamch=slamch\ -Ddlange=slange -Ddscal=sscal -Ddsyev=ssyev\ -Ddsyevd=ssyevd -Ddsyevx=ssyevx -Ddsymm=ssymm\ - -Ddsymv=ssymv -Ddsyr=ssyr -Ddsyrk=ssyrk -Dzcopy=ccopy\ - -Dzgemm=cgemm -Dzgemv=cgemv -Dzgerc=cgerc\ - -Dzgeru=cgeru -Dzscal=cscal -Dzsymm=csymm -Dzsyrk=csyrk -FCFLAGS = -D__T3E -D__FFTSG -D__parallel\ - -Ddcopy=scopy -Ddgbsv=sgbsv -Ddgecon=sgecon -Ddgemm=sgemm\ - -Ddgemv=sgemv -Ddger=sger -Ddgerfs=sgerfs -Ddgetrf=sgetrf\ - -Ddgetri=sgetri -Ddgetrs=sgetrs -Ddlamch=slamch\ - -Ddlange=slange -Ddscal=sscal -Ddsyev=ssyev\ - -Ddsyevd=ssyevd -Ddsyevx=ssyevx -Ddsymm=ssymm\ - -Ddsymv=ssymv -Ddsyr=ssyr -Ddsyrk=ssyrk -Dzcopy=ccopy\ - -Dzgemm=cgemm -Dzgemv=cgemv -Dzgerc=cgerc\ + -Ddsymv=ssymv -Ddsyr=ssyr -Ddsyrk=ssyrk\ + -Dzcopy=ccopy -Dzgemm=cgemm -Dzgemv=cgemv -Dzgerc=cgerc\ -Dzgeru=cgeru -Dzscal=cscal -Dzsymm=csymm -Dzsyrk=csyrk\ - -F -Racps -Xm -eIin -g -m2 -LDFLAGS = $(FCFLAGS) -L/u/krack/lib -LIBS = -llapack-dbg + -Ddgerv2d=sgerv2d -Ddgesd2d=sgesd2d\ + -Dpdgemm=psgemm -Dpdlamch=pslamch -Dpdsyevx=pssyevx\ + -Dpdsymm=pssymm -Dpdsyrk=pssyrk -Dpdtran=pstran +CPPFLAGS = -C $(DFLAGS) -P +FCFLAGS = $(DFLAGS) -F -Racps -Xm -eIin -g -m2 +LDFLAGS = $(FCFLAGS) +LIBS = OBJECTS_ARCHITECTURE = machine_t3e.o diff --git a/arch/CRAY-T3E.popt b/arch/CRAY-T3E.popt index 6289542fc3..29c74f11d8 100644 --- a/arch/CRAY-T3E.popt +++ b/arch/CRAY-T3E.popt @@ -5,26 +5,21 @@ FC = f90 -f free FC_fixed = f90 -f fixed LD = f90 AR = ar -r -CPPFLAGS = -C -D__T3E -D__FFTSG -D__parallel -P\ +DFLAGS = -D__T3E -D__FFTSG -D__parallel\ -Ddcopy=scopy -Ddgbsv=sgbsv -Ddgecon=sgecon -Ddgemm=sgemm\ -Ddgemv=sgemv -Ddger=sger -Ddgerfs=sgerfs -Ddgetrf=sgetrf\ -Ddgetri=sgetri -Ddgetrs=sgetrs -Ddlamch=slamch\ -Ddlange=slange -Ddscal=sscal -Ddsyev=ssyev\ -Ddsyevd=ssyevd -Ddsyevx=ssyevx -Ddsymm=ssymm\ - -Ddsymv=ssymv -Ddsyr=ssyr -Ddsyrk=ssyrk -Dzcopy=ccopy\ - -Dzgemm=cgemm -Dzgemv=cgemv -Dzgerc=cgerc\ - -Dzgeru=cgeru -Dzscal=cscal -Dzsymm=csymm -Dzsyrk=csyrk -FCFLAGS = -D__T3E -D__FFTSG -D__parallel\ - -Ddcopy=scopy -Ddgbsv=sgbsv -Ddgecon=sgecon -Ddgemm=sgemm\ - -Ddgemv=sgemv -Ddger=sger -Ddgerfs=sgerfs -Ddgetrf=sgetrf\ - -Ddgetri=sgetri -Ddgetrs=sgetrs -Ddlamch=slamch\ - -Ddlange=slange -Ddscal=sscal -Ddsyev=ssyev\ - -Ddsyevd=ssyevd -Ddsyevx=ssyevx -Ddsymm=ssymm\ - -Ddsymv=ssymv -Ddsyr=ssyr -Ddsyrk=ssyrk -Dzcopy=ccopy\ - -Dzgemm=cgemm -Dzgemv=cgemv -Dzgerc=cgerc\ + -Ddsymv=ssymv -Ddsyr=ssyr -Ddsyrk=ssyrk\ + -Dzcopy=ccopy -Dzgemm=cgemm -Dzgemv=cgemv -Dzgerc=cgerc\ -Dzgeru=cgeru -Dzscal=cscal -Dzsymm=csymm -Dzsyrk=csyrk\ - -F -O3 -Xm -LDFLAGS = $(FCFLAGS) -L/u/krack/lib -LIBS = -llapack-opt + -Ddgerv2d=sgerv2d -Ddgesd2d=sgesd2d\ + -Dpdgemm=psgemm -Dpdlamch=pslamch -Dpdsyevx=pssyevx\ + -Dpdsymm=pssymm -Dpdsyrk=pssyrk -Dpdtran=pstran +CPPFLAGS = -C $(DFLAGS) -P +FCFLAGS = $(DFLAGS) -F -O3 -Xm +LDFLAGS = $(FCFLAGS) +LIBS = OBJECTS_ARCHITECTURE = machine_t3e.o diff --git a/arch/CRAY-T3E.sdbg b/arch/CRAY-T3E.sdbg index a08b7c78a5..20e8aeaca1 100644 --- a/arch/CRAY-T3E.sdbg +++ b/arch/CRAY-T3E.sdbg @@ -5,26 +5,18 @@ FC = f90 -f free FC_fixed = f90 -f fixed LD = f90 AR = ar -r -CPPFLAGS = -C -D__T3E -D__FFTSG -P\ +DFLAGS = -D__T3E -D__FFTSG\ -Ddcopy=scopy -Ddgbsv=sgbsv -Ddgecon=sgecon -Ddgemm=sgemm\ -Ddgemv=sgemv -Ddger=sger -Ddgerfs=sgerfs -Ddgetrf=sgetrf\ -Ddgetri=sgetri -Ddgetrs=sgetrs -Ddlamch=slamch\ -Ddlange=slange -Ddscal=sscal -Ddsyev=ssyev\ -Ddsyevd=ssyevd -Ddsyevx=ssyevx -Ddsymm=ssymm\ - -Ddsymv=ssymv -Ddsyr=ssyr -Ddsyrk=ssyrk -Dzcopy=ccopy\ - -Dzgemm=cgemm -Dzgemv=cgemv -Dzgerc=cgerc\ + -Ddsymv=ssymv -Ddsyr=ssyr -Ddsyrk=ssyrk\ + -Dzcopy=ccopy -Dzgemm=cgemm -Dzgemv=cgemv -Dzgerc=cgerc\ -Dzgeru=cgeru -Dzscal=cscal -Dzsymm=csymm -Dzsyrk=csyrk -FCFLAGS = -D__T3E -D__FFTSG\ - -Ddcopy=scopy -Ddgbsv=sgbsv -Ddgecon=sgecon -Ddgemm=sgemm\ - -Ddgemv=sgemv -Ddger=sger -Ddgerfs=sgerfs -Ddgetrf=sgetrf\ - -Ddgetri=sgetri -Ddgetrs=sgetrs -Ddlamch=slamch\ - -Ddlange=slange -Ddscal=sscal -Ddsyev=ssyev\ - -Ddsyevd=ssyevd -Ddsyevx=ssyevx -Ddsymm=ssymm\ - -Ddsymv=ssymv -Ddsyr=ssyr -Ddsyrk=ssyrk -Dzcopy=ccopy\ - -Dzgemm=cgemm -Dzgemv=cgemv -Dzgerc=cgerc\ - -Dzgeru=cgeru -Dzscal=cscal -Dzsymm=csymm -Dzsyrk=csyrk\ - -F -Racps -eIin -g -m2 -LDFLAGS = $(FCFLAGS) -L/u/krack/lib -LIBS = -llapack-dbg +CPPFLAGS = -C $(DFLAGS) -P +FCFLAGS = $(DFLAGS) -F -Racps -eIin -g -m2 +LDFLAGS = $(FCFLAGS) +LIBS = OBJECTS_ARCHITECTURE = machine_t3e.o diff --git a/arch/CRAY-T3E.sopt b/arch/CRAY-T3E.sopt index b6c2c1bd7c..0293c6f46f 100644 --- a/arch/CRAY-T3E.sopt +++ b/arch/CRAY-T3E.sopt @@ -5,26 +5,18 @@ FC = f90 -f free FC_fixed = f90 -f fixed LD = f90 AR = ar -r -CPPFLAGS = -C -D__T3E -D__FFTSG -P\ +DFLAGS = -D__T3E -D__FFTSG\ -Ddcopy=scopy -Ddgbsv=sgbsv -Ddgecon=sgecon -Ddgemm=sgemm\ -Ddgemv=sgemv -Ddger=sger -Ddgerfs=sgerfs -Ddgetrf=sgetrf\ -Ddgetri=sgetri -Ddgetrs=sgetrs -Ddlamch=slamch\ -Ddlange=slange -Ddscal=sscal -Ddsyev=ssyev\ -Ddsyevd=ssyevd -Ddsyevx=ssyevx -Ddsymm=ssymm\ - -Ddsymv=ssymv -Ddsyr=ssyr -Ddsyrk=ssyrk -Dzcopy=ccopy\ - -Dzgemm=cgemm -Dzgemv=cgemv -Dzgerc=cgerc\ + -Ddsymv=ssymv -Ddsyr=ssyr -Ddsyrk=ssyrk\ + -Dzcopy=ccopy -Dzgemm=cgemm -Dzgemv=cgemv -Dzgerc=cgerc\ -Dzgeru=cgeru -Dzscal=cscal -Dzsymm=csymm -Dzsyrk=csyrk -FCFLAGS = -D__T3E -D__FFTSG\ - -Ddcopy=scopy -Ddgbsv=sgbsv -Ddgecon=sgecon -Ddgemm=sgemm\ - -Ddgemv=sgemv -Ddger=sger -Ddgerfs=sgerfs -Ddgetrf=sgetrf\ - -Ddgetri=sgetri -Ddgetrs=sgetrs -Ddlamch=slamch\ - -Ddlange=slange -Ddscal=sscal -Ddsyev=ssyev\ - -Ddsyevd=ssyevd -Ddsyevx=ssyevx -Ddsymm=ssymm\ - -Ddsymv=ssymv -Ddsyr=ssyr -Ddsyrk=ssyrk -Dzcopy=ccopy\ - -Dzgemm=cgemm -Dzgemv=cgemv -Dzgerc=cgerc\ - -Dzgeru=cgeru -Dzscal=cscal -Dzsymm=csymm -Dzsyrk=csyrk\ - -F -O3 -LDFLAGS = $(FCFLAGS) -L/u/krack/lib -LIBS = -llapack-opt +CPPFLAGS = -C $(DFLAGS) -P +FCFLAGS = $(DFLAGS) -F -O3 +LDFLAGS = $(FCFLAGS) +LIBS = OBJECTS_ARCHITECTURE = machine_t3e.o diff --git a/arch/Linux-i686-pgi.pdbg b/arch/Linux-i686-pgi.pdbg new file mode 100644 index 0000000000..cea8d77e4a --- /dev/null +++ b/arch/Linux-i686-pgi.pdbg @@ -0,0 +1,15 @@ +PERL = perl +CC = cc +CPP = cpp +FC = mpif90 -Mfree +FC_fixed = mpif90 -Mfixed +LD = mpif90 +AR = ar -r +DFLAGS = -D__PGI -D__FFTSG -D__parallel +CPPFLAGS = -C $(DFLAGS) -P +FCFLAGS = $(DFLAGS) -Mbounds -g +LDFLAGS = $(FCFLAGS) +LIBS = -lscalapack -lpblas -ltools -lblacsF77init -lblacs -lblacsF77init\ + -llapack -lf77blas -latlas + +OBJECTS_ARCHITECTURE = machine_pgi.o diff --git a/arch/Linux-i686-pgi.popt b/arch/Linux-i686-pgi.popt new file mode 100644 index 0000000000..6b089bc662 --- /dev/null +++ b/arch/Linux-i686-pgi.popt @@ -0,0 +1,15 @@ +PERL = perl +CC = cc +CPP = cpp +FC = mpif90 -Mfree +FC_fixed = mpif90 -Mfixed +LD = mpif90 +AR = ar -r +DFLAGS = -D__PGI -D__FFTSG -D__parallel +CPPFLAGS = -C $(DFLAGS) -P +FCFLAGS = $(DFLAGS) -fast +LDFLAGS = $(FCFLAGS) +LIBS = -lscalapack -lpblas -ltools -lblacsF77init -lblacs -lblacsF77init\ + -llapack -lf77blas -latlas + +OBJECTS_ARCHITECTURE = machine_pgi.o diff --git a/arch/Linux-i686-pgi.sdbg b/arch/Linux-i686-pgi.sdbg index 06811aba01..fc674602bb 100644 --- a/arch/Linux-i686-pgi.sdbg +++ b/arch/Linux-i686-pgi.sdbg @@ -5,8 +5,9 @@ FC = pgf90 -Mfree FC_fixed = pgf90 -Mfixed LD = pgf90 AR = ar -r -CPPFLAGS = -C -D__PGI -D__FFTSG -P -FCFLAGS = -Mbounds -g -D__PGI -D__FFTSG +DFLAGS = -D__PGI -D__FFTSG +CPPFLAGS = -C $(DFLAGS) -P +FCFLAGS = $(DFLAGS) -Mbounds -g LDFLAGS = $(FCFLAGS) LIBS = -llapack -lblas diff --git a/arch/Linux-i686-pgi.sopt b/arch/Linux-i686-pgi.sopt index 9ea85f902d..83fcea7b82 100644 --- a/arch/Linux-i686-pgi.sopt +++ b/arch/Linux-i686-pgi.sopt @@ -5,8 +5,9 @@ FC = pgf90 -Mfree FC_fixed = pgf90 -Mfixed LD = pgf90 AR = ar -r -CPPFLAGS = -C -D__PGI -D__FFTSG -P -FCFLAGS = -fast -D__PGI -D__FFTSG +DFLAGS = -D__PGI -D__FFTSG +CPPFLAGS = -C $(DFLAGS) -P +FCFLAGS = $(DFLAGS) -fast LDFLAGS = $(FCFLAGS) LIBS = -llapack -lblas diff --git a/src/OBJECTDEFS b/src/OBJECTDEFS index d13c8d304b..881f7683bb 100644 --- a/src/OBJECTDEFS +++ b/src/OBJECTDEFS @@ -5,6 +5,7 @@ OBJECTS_GENERIC =\ ai_nuclear.o\ ai_overlap.o\ ai_overlap_ppl.o\ + ai_overlap3.o\ ai_verfc.o\ amoeba.o\ ao_types.o\ @@ -15,6 +16,7 @@ OBJECTS_GENERIC =\ band.o\ basis_set_types.o\ bessel_lib.o\ + blacs.o\ brillouin.o\ cell_parameters.o\ cntl_input.o\ diff --git a/src/blacs.F b/src/blacs.F new file mode 100644 index 0000000000..b7746104b4 --- /dev/null +++ b/src/blacs.F @@ -0,0 +1,1958 @@ +!-----------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright (C) 1999 MPI fuer Festkoerperforschung, Stuttgart ! +!-----------------------------------------------------------------------------! +!!****** cp2k/blacs [1.0] * +!! +!! NAME +!! blacs +!! +!! FUNCTION +!! BLACS +!! +!! AUTHOR +!! Matthias Krack (22.05.2001) +!! +!! MODIFICATION HISTORY +!! none +!! +!! SOURCE +!****************************************************************************** + +MODULE blacs + +! ***************************************************************************** + + USE kinds, ONLY: wp => dp + + USE global_types, ONLY: global_environment_type + USE mathlib, ONLY: symmetrize_matrix + USE matrix_types, ONLY: first_block_node,& + get_block_node,& + get_matrix_info,& + next_block_node,& + real_block_node_type,& + real_matrix_type + USE memory_utilities, ONLY: reallocate + USE message_passing, ONLY: mp_bcast,mp_max,mp_sum + USE termination, ONLY: stop_program + USE string_utilities, ONLY: compress + USE timings, ONLY: timeset,timestop + + IMPLICIT NONE + + PRIVATE + + TYPE blacs_matrix_block_type + PRIVATE + INTEGER :: ncol_local,nrow_local + REAL(wp), DIMENSION(:,:), POINTER :: block + END TYPE blacs_matrix_block_type + + TYPE blacs_matrix_type + PRIVATE + CHARACTER(LEN=60) :: name + INTEGER :: context,& + ncol_block,& + ncol_global,& + nrow_block,& + nrow_global + INTEGER, DIMENSION(:), POINTER :: descriptor + TYPE(blacs_matrix_block_type), DIMENSION(:,:), POINTER :: p + END TYPE blacs_matrix_type + +! *** Public data types *** + + PUBLIC :: blacs_matrix_type + +! *** Public subroutines *** + + PUBLIC :: allocate_blacs_matrix,& + blacs_add,& + blacs_gemm,& + blacs_maxval,& + blacs_set_all,& + blacs_set_element,& + blacs_syevx,& + blacs_symm,& + blacs_syrk,& + blacs_trace,& + copy_blacs_to_blacs_matrix,& + copy_blacs_to_full_matrix,& + copy_blacs_to_sparse_matrix,& + copy_sparse_to_blacs_matrix,& + deallocate_blacs_matrix,& + finish_blacs,& + get_blacs_info,& + get_blacs_matrix_info,& + power_blacs_matrix,& + read_blacs_matrix,& + replicate_blacs_matrix,& + start_blacs,& + symmetrise_blacs_matrix,& + write_blacs_matrix + +! ***************************************************************************** + +CONTAINS + +! ***************************************************************************** + + SUBROUTINE allocate_blacs_matrix(new_matrix,nrow_global,ncol_global,& + nrow_block,ncol_block,name,context,globenv) + +! Purpose: Allocate a new distributed BLACS matrix. + +! History: - Creation (23.05.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(OUT) :: new_matrix + TYPE(global_environment_type), INTENT(IN) :: globenv + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER, INTENT(IN) :: context,ncol_block,& + ncol_global,nrow_block,& + nrow_global + +! *** Local parameters *** + + CHARACTER(LEN=*), PARAMETER :: routine =& + "SUBROUTINE allocate_blacs_matrix (MODULE blacs)" + +! *** Local variables *** + + CHARACTER(LEN=40) :: message + INTEGER :: group,ierror,ipcol,ipe,iprow,mype,mypcol,myprow,& + ncol_local,npcol,npe,nprow,nrow_local,output_unit,& + source + LOGICAL :: ionode + + INTEGER, DIMENSION(:), POINTER :: pcol,prow + +#if defined(__parallel) + INTEGER, EXTERNAL :: blacs_pnum,numroc + +#endif +! --------------------------------------------------------------------------- + + group = globenv%group + ionode = globenv%ionode + output_unit = globenv%scr + source = globenv%source + + new_matrix%name = name +#if defined(__parallel) + + CALL blacs_pinfo(mype,npe) + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + + new_matrix%context = context + new_matrix%nrow_global = nrow_global + new_matrix%ncol_global = ncol_global + + new_matrix%nrow_block = MIN(nrow_block,nrow_global/nprow,& + ncol_block,ncol_global/npcol) + new_matrix%ncol_block = new_matrix%nrow_block + + IF ((new_matrix%nrow_block == 0).OR.& + (new_matrix%ncol_block == 0)) THEN + CALL stop_program(routine,"More processes than matrix elements",globenv) + END IF + + nrow_local = numroc(nrow_global,new_matrix%nrow_block,myprow,source,nprow) + ncol_local = numroc(ncol_global,new_matrix%ncol_block,mypcol,source,npcol) + + NULLIFY (prow,pcol) + prow => reallocate(prow,0,npe-1) + pcol => reallocate(pcol,0,npe-1) + + prow(mype) = nrow_local + pcol(mype) = ncol_local + + CALL mp_sum(prow,group) + CALL mp_sum(pcol,group) + + NULLIFY (new_matrix%descriptor) + new_matrix%descriptor => reallocate(new_matrix%descriptor,1,9) + + IF (ionode) THEN + WRITE (UNIT=output_unit,FMT="(/,/,T2,A,/,/,T3,A,/,/,(T3,A,I6))")& + "BLACS INFORMATION (BLACS matrix allocation)",& + "Matrix name: "//TRIM(name),& + "Number of rows of the global matrix: ",new_matrix%nrow_global,& + "Number of columns of the global matrix: ",new_matrix%ncol_global,& + "Number of rows of a matrix block: ",new_matrix%nrow_block,& + "Number of columns of a matrix block: ",new_matrix%ncol_block + WRITE (UNIT=output_unit,FMT="(/,T4,A,/)")& + "PE block rows block columns rows columns" + WRITE (UNIT=output_unit,FMT="(I5,T16,I6,T32,I6,T42,I6,T52,I6)")& + (ipe,prow(ipe)/new_matrix%nrow_block,pcol(ipe)/new_matrix%nrow_block,& + prow(ipe),pcol(ipe),ipe=0,npe-1) + END IF + + ALLOCATE (new_matrix%p(0:nprow-1,0:npcol-1)) + + DO iprow=0,nprow-1 + DO ipcol=0,npcol-1 + NULLIFY (new_matrix%p(iprow,ipcol)%block) + ipe = blacs_pnum(context,iprow,ipcol) + new_matrix%p(iprow,ipcol)%nrow_local = prow(ipe) + new_matrix%p(iprow,ipcol)%ncol_local = pcol(ipe) + END DO + END DO + + new_matrix%p(myprow,mypcol)%block =>& + reallocate(new_matrix%p(myprow,mypcol)%block,& + 1,nrow_local,& + 1,ncol_local) + + CALL descinit(new_matrix%descriptor,new_matrix%nrow_global,& + new_matrix%ncol_global,new_matrix%nrow_block,& + new_matrix%ncol_block,source,source,context,nrow_local,& + ierror) + + IF (ierror /= 0) THEN + WRITE (UNIT=message,FMT="(A,I6)") "Error in descinit: ierror = ",ierror + CALL compress(message) + CALL stop_program(routine,message,globenv) + END IF + + DEALLOCATE (prow,pcol) + +#else + + new_matrix%context = 0 + new_matrix%nrow_block = nrow_global + new_matrix%ncol_block = ncol_global + new_matrix%nrow_global = nrow_global + new_matrix%ncol_global = ncol_global + NULLIFY (new_matrix%descriptor) + new_matrix%descriptor => reallocate(new_matrix%descriptor,1,9) + ALLOCATE (new_matrix%p(source:source,source:source)) + NULLIFY (new_matrix%p(source,source)%block) + new_matrix%p(source,source)%nrow_local = nrow_global + new_matrix%p(source,source)%ncol_local = ncol_global + new_matrix%p(source,source)%block =>& + reallocate(new_matrix%p(source,source)%block,& + 1,new_matrix%p(source,source)%nrow_local,& + 1,new_matrix%p(source,source)%ncol_local) + +#endif + END SUBROUTINE allocate_blacs_matrix + +! ***************************************************************************** + + SUBROUTINE blacs_add(alpha,matrix_a,beta,matrix_b,context,globenv) + +! Purpose: Scale and add two BLACS matrices (a <- alpha*a + beta*b). + +! History: - Creation (11.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(INOUT) :: matrix_a + TYPE(blacs_matrix_type), INTENT(IN) :: matrix_b + TYPE(global_environment_type), INTENT(IN) :: globenv + REAL(wp), INTENT(IN) :: alpha,beta + INTEGER, INTENT(IN) :: context + +! *** Local variables *** + + INTEGER :: handle,mypcol,myprow,npcol,nprow,source + + REAL(wp), DIMENSION(:,:), POINTER :: a,b + +! --------------------------------------------------------------------------- + + CALL timeset("blacs_add","I","",handle) + + source = globenv%source +#if defined(__parallel) + + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + +#else + + myprow = source + mypcol = source + +#endif + a => matrix_a%p(myprow,mypcol)%block + b => matrix_b%p(myprow,mypcol)%block + + IF (alpha == 0.0_wp) THEN + IF (beta == 0.0_wp) THEN + a(:,:) = 0.0_wp + ELSE IF (beta == 1.0_wp) THEN + a(:,:) = b(:,:) + ELSE + a(:,:) = beta*b(:,:) + END IF + ELSE IF (beta == 0.0_wp) THEN + IF (alpha == 1.0_wp) THEN + RETURN + ELSE + a(:,:) = alpha*a(:,:) + END IF + ELSE IF (alpha == 1.0_wp) THEN + IF (beta == 1.0_wp) THEN + a(:,:) = a(:,:) + b(:,:) + ELSE + a(:,:) = a(:,:) + beta*b(:,:) + END IF + ELSE IF (beta == 1.0_wp) THEN + a(:,:) = alpha*a(:,:) + b(:,:) + ELSE + a(:,:) = alpha*a(:,:) + beta*b(:,:) + END IF + + CALL timestop(0.0_wp,handle) + + END SUBROUTINE blacs_add + +! ***************************************************************************** + + SUBROUTINE blacs_gemm(transa,transb,m,n,k,alpha,matrix_a,matrix_b,beta,& + matrix_c,context,globenv) + +! Purpose: BLACS interface to the BLAS routine dgemm. + +! History: - Creation (07.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(IN) :: matrix_a,matrix_b + TYPE(blacs_matrix_type), INTENT(INOUT) :: matrix_c + TYPE(global_environment_type), INTENT(IN) :: globenv + CHARACTER(LEN=1), INTENT(IN) :: transa,transb + REAL(wp), INTENT(IN) :: alpha,beta + INTEGER, INTENT(IN) :: context,k,m,n + +! *** Local variables *** + + INTEGER :: handle,lda,ldb,ldc,mypcol,myprow,npcol,nprow,source + + INTEGER, DIMENSION(:), POINTER :: desca,descb,descc + REAL(wp), DIMENSION(:,:), POINTER :: a,b,c + +! --------------------------------------------------------------------------- + + CALL timeset("blacs_gemm","I","",handle) + + source = globenv%source +#if defined(__parallel) + + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + + a => matrix_a%p(myprow,mypcol)%block + desca => matrix_a%descriptor + b => matrix_b%p(myprow,mypcol)%block + descb => matrix_b%descriptor + c => matrix_c%p(myprow,mypcol)%block + descc => matrix_c%descriptor + + CALL pdgemm(transa,transb,m,n,k,alpha,a,1,1,desca,b,1,1,descb,beta,c,1,1,& + descc) + +#else + + a => matrix_a%p(source,source)%block + b => matrix_b%p(source,source)%block + c => matrix_c%p(source,source)%block + + lda = matrix_a%nrow_global + ldb = matrix_b%nrow_global + ldc = matrix_c%nrow_global + + CALL dgemm(transa,transb,m,n,k,alpha,a,lda,b,ldb,beta,c,ldc) + +#endif + CALL timestop(0.0_wp,handle) + + END SUBROUTINE blacs_gemm + +! ***************************************************************************** + + SUBROUTINE blacs_maxval(matrix,a_max,context,globenv) + +! Purpose: Get the maximum absolute element of a BLACS matrix. + +! History: - Creation (11.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(IN) :: matrix + TYPE(global_environment_type), INTENT(IN) :: globenv + REAL(wp), INTENT(OUT) :: a_max + INTEGER, INTENT(IN) :: context + +! *** Local variables *** + + INTEGER :: group,handle,mypcol,myprow,npcol,nprow,source + +! --------------------------------------------------------------------------- + + CALL timeset("blacs_maxval","I","",handle) + + source = globenv%source + group = globenv%group +#if defined(__parallel) + + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + +#else + + myprow = source + mypcol = source + +#endif + a_max = MAXVAL(ABS(matrix%p(myprow,mypcol)%block)) + + CALL mp_max(a_max,group) + + CALL timestop(0.0_wp,handle) + + END SUBROUTINE blacs_maxval + +! ***************************************************************************** + + SUBROUTINE blacs_set_all(matrix,alpha,context,globenv) + +! Purpose: Set the BLACS matrix elements to alpha. + +! History: - Creation (12.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(INOUT) :: matrix + TYPE(global_environment_type), INTENT(IN) :: globenv + REAL(wp), INTENT(IN) :: alpha + INTEGER, INTENT(IN) :: context + +! *** Local variables *** + + INTEGER :: mypcol,myprow,npcol,nprow,source + +! --------------------------------------------------------------------------- + + source = globenv%source +#if defined(__parallel) + + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + + matrix%p(myprow,mypcol)%block(:,:) = alpha + +#else + + matrix%p(source,source)%block(:,:) = alpha + +#endif + END SUBROUTINE blacs_set_all + +! ***************************************************************************** + + SUBROUTINE blacs_set_element(matrix,irow_global,icol_global,alpha,context,& + globenv) + +! Purpose: Set the BLACS matrix element (irow_global,icol_global) to alpha. + +! History: - Creation (08.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(INOUT) :: matrix + TYPE(global_environment_type), INTENT(IN) :: globenv + REAL(wp), INTENT(IN) :: alpha + INTEGER, INTENT(IN) :: context,icol_global,& + irow_global + +! *** Local variables *** + + INTEGER :: icol_local,ipcol,iprow,irow_local,mypcol,myprow,npcol,nprow,& + source + + INTEGER, DIMENSION(:), POINTER :: desca + REAL(wp), DIMENSION(:,:), POINTER :: a + +! --------------------------------------------------------------------------- + + source = globenv%source +#if defined(__parallel) + + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + + a => matrix%p(myprow,mypcol)%block + desca => matrix%descriptor + + CALL infog2l(irow_global,icol_global,desca,nprow,npcol,myprow,mypcol,& + irow_local,icol_local,iprow,ipcol) + + IF ((iprow == myprow).AND.(ipcol == mypcol)) THEN + a(irow_local,icol_local) = alpha + END IF + +#else + + matrix%p(source,source)%block(irow_global,icol_global) = alpha + +#endif + END SUBROUTINE blacs_set_element + +! ***************************************************************************** + + SUBROUTINE blacs_syevx(matrix,eigenvectors,eigenvalues,neig,work_syevx,& + context,globenv) + +! Purpose: Diagonalise the symmetric n by n matrix using the LAPACK library. + +! History: - Creation (06.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(INOUT) :: matrix + TYPE(blacs_matrix_type), INTENT(OUT) :: eigenvectors + TYPE(global_environment_type), INTENT(IN) :: globenv + REAL(wp), INTENT(IN) :: work_syevx + INTEGER, INTENT(IN) :: context,neig + REAL(wp), DIMENSION(:), INTENT(OUT) :: eigenvalues + +! *** Local parameters *** + + CHARACTER(LEN=*), PARAMETER :: routine =& + "SUBROUTINE blas_syevx (MODULE blacs)" + REAL(wp), PARAMETER :: orfac = -1.0_wp,& + vl = 0.0_wp,& + vu = 0.0_wp + +! *** Local variables *** + + REAL(wp) :: abstol + INTEGER :: handle,info,liwork,lwork,m,mypcol,myprow,n,nb,nn,np0,npcol,& + npe,nprow,nq0,nz,output_unit,source + LOGICAL :: ionode + + REAL(wp), DIMENSION(:), POINTER :: gap,work + INTEGER, DIMENSION(:), POINTER :: desca,descz,iclustr,ifail,iwork + REAL(wp), DIMENSION(:,:), POINTER :: a,z + +#if defined(__parallel) + REAL(wp), EXTERNAL :: pdlamch + INTEGER, EXTERNAL :: iceil,numroc +#else + REAL(wp), EXTERNAL :: dlamch + INTEGER, EXTERNAL :: ilaenv +#endif +! --------------------------------------------------------------------------- + + CALL timeset("blacs_syevx","I","",handle) + + ionode = globenv%ionode + output_unit = globenv%scr + source = globenv%source + + n = matrix%nrow_global +#if defined(__parallel) + + IF (matrix%nrow_block /= matrix%ncol_block) THEN + CALL stop_program(routine,"Invalid blocksize (no square blocks)",globenv) + END IF + + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + + a => matrix%p(myprow,mypcol)%block + desca => matrix%descriptor + z => eigenvectors%p(myprow,mypcol)%block + descz => eigenvectors%descriptor + +! *** Get the optimal work storage size *** + + npe = nprow*npcol + nb = matrix%nrow_block + nn = MAX(n,nb,2) + np0 = numroc(nn,nb,0,0,nprow) + nq0 = MAX(numroc(nn,nb,0,0,npcol),nb) + + lwork = 5*n + MAX(5*nn,np0*nq0) + iceil(neig,npe)*nn + 2*nb*nb +& + INT(work_syevx*REAL((neig - 1)*n,wp)) + liwork = MAX(3*n + npe + 1,4*n,14) + 2*n + + NULLIFY (gap,iclustr,ifail,iwork,work) + gap => reallocate(gap,1,npe) + iclustr => reallocate(iclustr,1,2*npe) + ifail => reallocate(ifail,1,n) + iwork => reallocate(iwork,1,liwork) + work => reallocate(work,1,lwork) + +! *** Set the absolute error tolerance for the eigenvalues *** + + abstol = pdlamch(context,"U") + +! *** Diagonalise matrix *** + + CALL pdsyevx("V","I","U",n,a,1,1,desca,vl,vu,1,neig,abstol,m,nz,& + eigenvalues,orfac,z,1,1,descz,work,lwork,iwork,liwork,ifail,& + iclustr,gap,info) + +! *** Error handling *** + + IF (info /= 0) THEN + IF (ionode) THEN + WRITE (unit=output_unit,FMT="(/,(T3,A,T12,1X,I10))")& + "info = ",info,& + "lwork = ",lwork,& + "liwork = ",liwork,& + "nz = ",nz + IF (info > 0) THEN + WRITE (unit=output_unit,FMT="(/,T3,A,(T12,6(1X,I10)))")& + "ifail = ",ifail + WRITE (unit=output_unit,FMT="(/,T3,A,(T12,6(1X,I10)))")& + "iclustr = ",iclustr + WRITE (unit=output_unit,FMT="(/,T3,A,(T12,6(1X,E10.3)))")& + "gap = ",gap + END IF + END IF + CALL stop_program(routine,"Error in pdsyevx",globenv) + END IF + +! *** Release work storage *** + + DEALLOCATE (gap,iclustr,ifail,iwork,work) + +#else + + a => matrix%p(source,source)%block + z => eigenvectors%p(source,source)%block + +! *** Get the optimal work storage size *** + + nb = MAX(ilaenv(1,"DSYTRD","U",n,-1,-1,-1),& + ilaenv(1,"DORMTR","U",n,-1,-1,-1),8*n) + + lwork = (nb + 3)*n + liwork = 5*n + + NULLIFY (ifail,iwork,work) + ifail => reallocate(ifail,1,n) + iwork => reallocate(iwork,1,liwork) + work => reallocate(work,1,lwork) + +! *** Set the absolute error tolerance for the eigenvalues *** + + abstol = 2.0_wp*dlamch("S") + +! *** Diagonalise matrix *** + + CALL dsyevx("V","I","U",n,a,n,vl,vu,1,neig,abstol,m,eigenvalues,z,n,work,& + lwork,iwork,ifail,info) + +! *** Error handling *** + + IF (info /= 0) CALL stop_program(routine,"Error in dsyevx",globenv) + +! *** Release work storage *** + + DEALLOCATE (ifail,iwork,work) + +#endif + CALL timestop(0.0_wp,handle) + + END SUBROUTINE blacs_syevx + +! ***************************************************************************** + + SUBROUTINE blacs_symm(side,uplo,m,n,alpha,matrix_a,matrix_b,beta,matrix_c,& + context,globenv) + +! Purpose: BLACS interface to the BLAS routine dsymm. + +! History: - Creation (07.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(IN) :: matrix_a,matrix_b + TYPE(blacs_matrix_type), INTENT(INOUT) :: matrix_c + TYPE(global_environment_type), INTENT(IN) :: globenv + CHARACTER(LEN=1), INTENT(IN) :: side,uplo + REAL(wp), INTENT(IN) :: alpha,beta + INTEGER, INTENT(IN) :: context,m,n + +! *** Local variables *** + + INTEGER :: handle,lda,ldb,ldc,mypcol,myprow,npcol,nprow,source + + INTEGER, DIMENSION(:), POINTER :: desca,descb,descc + REAL(wp), DIMENSION(:,:), POINTER :: a,b,c + +! --------------------------------------------------------------------------- + + CALL timeset("blacs_symm","I","",handle) + + source = globenv%source +#if defined(__parallel) + + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + + a => matrix_a%p(myprow,mypcol)%block + desca => matrix_a%descriptor + b => matrix_b%p(myprow,mypcol)%block + descb => matrix_b%descriptor + c => matrix_c%p(myprow,mypcol)%block + descc => matrix_c%descriptor + + CALL pdsymm(side,uplo,m,n,alpha,a,1,1,desca,b,1,1,descb,beta,c,1,1,descc) + +#else + + a => matrix_a%p(source,source)%block + b => matrix_b%p(source,source)%block + c => matrix_c%p(source,source)%block + + lda = matrix_a%nrow_global + ldb = matrix_b%nrow_global + ldc = matrix_c%nrow_global + + CALL dsymm(side,uplo,m,n,alpha,a,lda,b,ldb,beta,c,ldc) + +#endif + CALL timestop(0.0_wp,handle) + + END SUBROUTINE blacs_symm + +! ***************************************************************************** + + SUBROUTINE blacs_syrk(uplo,trans,k,alpha,matrix_a,beta,matrix_c,context,& + globenv) + +! Purpose: BLACS interface to the BLAS routine dsyrk. + +! History: - Creation (07.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(IN) :: matrix_a + TYPE(blacs_matrix_type), INTENT(INOUT) :: matrix_c + TYPE(global_environment_type), INTENT(IN) :: globenv + CHARACTER(LEN=1), INTENT(IN) :: trans,uplo + REAL(wp), INTENT(IN) :: alpha,beta + INTEGER, INTENT(IN) :: context,k + +! *** Local variables *** + + INTEGER :: handle,lda,ldc,mypcol,myprow,n,npcol,nprow,source + + INTEGER, DIMENSION(:), POINTER :: desca,descc + REAL(wp), DIMENSION(:,:), POINTER :: a,c + +! --------------------------------------------------------------------------- + + CALL timeset("blacs_syrk","I","",handle) + + source = globenv%source + n = matrix_a%nrow_global +#if defined(__parallel) + + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + + a => matrix_a%p(myprow,mypcol)%block + desca => matrix_a%descriptor + c => matrix_c%p(myprow,mypcol)%block + descc => matrix_c%descriptor + + CALL pdsyrk(uplo,trans,n,k,alpha,a,1,1,desca,beta,c,1,1,descc) + +#else + + a => matrix_a%p(source,source)%block + c => matrix_c%p(source,source)%block + + lda = matrix_a%nrow_global + ldc = matrix_c%nrow_global + + CALL dsyrk(uplo,trans,n,k,alpha,a,lda,beta,c,ldc) + +#endif + CALL timestop(0.0_wp,handle) + + END SUBROUTINE blacs_syrk + +! ***************************************************************************** + + SUBROUTINE blacs_trace(matrix_a,matrix_b,trace,context,globenv) + +! Purpose: Calculate the trace of the product of two BLACS matrices. + +! History: - Creation (11.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(IN) :: matrix_a,matrix_b + TYPE(global_environment_type), INTENT(IN) :: globenv + REAL(wp), INTENT(OUT) :: trace + INTEGER, INTENT(IN) :: context + +! *** Local variables *** + + INTEGER :: group,handle,icol_local,irow_local,mypcol,myprow,ncol_local,& + npcol,nprow,nrow_local,source + + REAL(wp), DIMENSION(:,:), POINTER :: a,b + +! --------------------------------------------------------------------------- + + CALL timeset("blacs_trace","I","",handle) + + group = globenv%group + source = globenv%source +#if defined(__parallel) + + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + +#else + + myprow = source + mypcol = source + +#endif + a => matrix_a%p(myprow,mypcol)%block + b => matrix_b%p(myprow,mypcol)%block + + nrow_local = matrix_a%p(myprow,mypcol)%nrow_local + ncol_local = matrix_b%p(myprow,mypcol)%ncol_local + + trace = 0.0_wp + + DO icol_local=1,ncol_local + DO irow_local=1,nrow_local + trace = trace + a(irow_local,icol_local)*b(irow_local,icol_local) + END DO + END DO + + CALL mp_sum(trace,group) + + CALL timestop(0.0_wp,handle) + + END SUBROUTINE blacs_trace + +! ***************************************************************************** + + SUBROUTINE copy_blacs_to_blacs_matrix(source_matrix,target_matrix) + +! Purpose: Copy BLACS matrix to a BLACS matrix of the same type. + +! History: - Creation (08.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(IN) :: source_matrix + TYPE(blacs_matrix_type), INTENT(OUT) :: target_matrix + +! *** Local variables *** + + INTEGER :: ipcol,iprow,npcol,nprow + +! --------------------------------------------------------------------------- + + nprow = SIZE(source_matrix%p,1) + npcol = SIZE(source_matrix%p,2) + + DO iprow=0,nprow-1 + DO ipcol=0,npcol-1 + IF (ASSOCIATED(source_matrix%p(iprow,ipcol)%block)) THEN + IF (.NOT.ASSOCIATED(target_matrix%p(iprow,ipcol)%block)) THEN + target_matrix%p(iprow,ipcol)%block =>& + reallocate(target_matrix%p(iprow,ipcol)%block,& + 1,target_matrix%p(iprow,ipcol)%nrow_local,& + 1,target_matrix%p(iprow,ipcol)%ncol_local) + END IF + target_matrix%p(iprow,ipcol)%block(:,:) =& + source_matrix%p(iprow,ipcol)%block(:,:) + ELSE + IF (ASSOCIATED(target_matrix%p(iprow,ipcol)%block)) THEN + DEALLOCATE (target_matrix%p(iprow,ipcol)%block) + END IF + END IF + END DO + END DO + + END SUBROUTINE copy_blacs_to_blacs_matrix + +! ***************************************************************************** + + SUBROUTINE copy_blacs_to_full_matrix(blacs_matrix,full_matrix,context,& + globenv) + +! Purpose: Copy a BLACS matrix to a full matrix. + +! History: - Creation (18.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(INOUT) :: blacs_matrix + TYPE(global_environment_type), INTENT(IN) :: globenv + INTEGER, INTENT(IN) :: context + REAL(wp), DIMENSION(:,:), POINTER :: full_matrix + +! *** Local variables *** + + INTEGER :: handle,icol_global,icol_local,ipcol,ipe,iprow,irow_global,& + irow_local,mypcol,mype,myprow,ncol_block,ncol_global,& + ncol_local,npcol,npe,nprow,nrow_block,nrow_global,nrow_local,& + source + LOGICAL :: ionode + + REAL(wp), DIMENSION(:,:), POINTER :: blacs_block + +#if defined(__parallel) + INTEGER, EXTERNAL :: blacs_pnum,indxl2g + +#endif +! --------------------------------------------------------------------------- + + CALL timeset("copy_blacs_to_full_matrix","I","",handle) + + ionode = globenv%ionode + source = globenv%source + + nrow_global = blacs_matrix%nrow_global + ncol_global = blacs_matrix%ncol_global + + IF (ionode) THEN + full_matrix => reallocate(full_matrix,1,nrow_global,1,ncol_global) + END IF +#if defined(__parallel) + + CALL blacs_pinfo(mype,npe) + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + + nrow_block = blacs_matrix%nrow_block + ncol_block = blacs_matrix%ncol_block + + DO iprow=0,nprow-1 + DO ipcol=0,npcol-1 + + ipe = blacs_pnum(context,iprow,ipcol) + + nrow_local = blacs_matrix%p(iprow,ipcol)%nrow_local + ncol_local = blacs_matrix%p(iprow,ipcol)%ncol_local + + IF (ionode) THEN + + IF (ipe /= mype) THEN + blacs_matrix%p(iprow,ipcol)%block =>& + reallocate(blacs_matrix%p(iprow,ipcol)%block,& + 1,blacs_matrix%p(iprow,ipcol)%nrow_local,& + 1,blacs_matrix%p(iprow,ipcol)%ncol_local) + CALL dgerv2d(context,nrow_local,ncol_local,& + blacs_matrix%p(iprow,ipcol)%block,nrow_local,& + iprow,ipcol) + END IF + + blacs_block => blacs_matrix%p(iprow,ipcol)%block + + DO icol_local=1,ncol_local + icol_global = indxl2g(icol_local,ncol_block,ipcol,source,npcol) + DO irow_local=1,nrow_local + irow_global = indxl2g(irow_local,nrow_block,iprow,source,nprow) + full_matrix(irow_global,icol_global) = blacs_block(irow_local,& + icol_local) + END DO + END DO + + IF (ipe /= mype) DEALLOCATE (blacs_matrix%p(iprow,ipcol)%block) + + ELSE + + IF (ipe == mype) THEN + CALL dgesd2d(context,nrow_local,ncol_local,& + blacs_matrix%p(iprow,ipcol)%block,nrow_local,& + source,source) + END IF + + END IF + + CALL blacs_barrier(context,"A") + + END DO + END DO + +#else + + full_matrix(:,:) = blacs_matrix%p(source,source)%block(:,:) + +#endif + CALL timestop(0.0_wp,handle) + + END SUBROUTINE copy_blacs_to_full_matrix + +! ***************************************************************************** + + SUBROUTINE copy_blacs_to_sparse_matrix(blacs_matrix,sparse_matrix,context,& + globenv) + +! Purpose: Copy a BLACS matrix to a sparse matrix. The BLACS matrix blocks +! are deallocated during the copy procedure. + +! History: - Creation (06.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(INOUT) :: blacs_matrix + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(real_matrix_type), POINTER :: sparse_matrix + INTEGER, INTENT(IN) :: context + +! *** Local variables *** + + TYPE(real_block_node_type), POINTER :: block_node + + INTEGER :: group,handle,iblock_col,iblock_row,icol,icol_global,icol_local,& + ipcol,ipe,iprow,irow,irow_global,irow_local,jpcol,jprow,mypcol,& + mype,myprow,nblock_row,ncol_block,npcol,npe,nprow,nrow_block,& + source + + INTEGER, DIMENSION(:), POINTER :: first_col,first_row,last_col,last_row + REAL(wp), DIMENSION(:,:), POINTER :: blacs_block,sparse_block + +#if defined(__parallel) + INTEGER, EXTERNAL :: blacs_pnum,indxg2l,indxg2p + +#endif +! --------------------------------------------------------------------------- + + CALL timeset("copy_blacs_to_sparse_matrix","I","",handle) + + group = globenv%group + source = globenv%source + + CALL get_matrix_info(matrix=sparse_matrix,& + nblock_row=nblock_row,& + first_row=first_row,& + first_col=first_col,& + last_row=last_row,& + last_col=last_col) +#if defined(__parallel) + + CALL blacs_pinfo(mype,npe) + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + + nrow_block = blacs_matrix%nrow_block + ncol_block = blacs_matrix%ncol_block + + DO iprow=0,nprow-1 + DO ipcol=0,npcol-1 + + ipe = blacs_pnum(context,iprow,ipcol) + + IF (ipe /= mype) THEN + blacs_matrix%p(iprow,ipcol)%block =>& + reallocate(blacs_matrix%p(iprow,ipcol)%block,& + 1,blacs_matrix%p(iprow,ipcol)%nrow_local,& + 1,blacs_matrix%p(iprow,ipcol)%ncol_local) + END IF + + blacs_block => blacs_matrix%p(iprow,ipcol)%block + + CALL mp_bcast(blacs_block,ipe,group) + + DO iblock_row=1,nblock_row + + block_node => first_block_node(matrix=sparse_matrix,& + block_row=iblock_row) + + DO WHILE (ASSOCIATED(block_node)) + + CALL get_block_node(block_node=block_node,& + block_col=iblock_col,& + block=sparse_block) + + icol = 1 + + DO icol_global=first_col(iblock_col),last_col(iblock_col) + + jpcol = indxg2p(icol_global,ncol_block,mypcol,source,npcol) + + IF (jpcol == ipcol) THEN + + icol_local = indxg2l(icol_global,ncol_block,mypcol,source,& + npcol) + + irow = 1 + + DO irow_global=first_row(iblock_row),last_row(iblock_row) + + jprow = indxg2p(irow_global,nrow_block,myprow,source,nprow) + + IF (jprow == iprow) THEN + + irow_local = indxg2l(irow_global,nrow_block,myprow,source,& + nprow) + + sparse_block(irow,icol) = blacs_block(irow_local,& + icol_local) + + END IF + + irow = irow + 1 + + END DO + + END IF + + icol = icol + 1 + + END DO + + block_node => next_block_node(block_node) + + END DO + + END DO + + IF (ipe /= mype) DEALLOCATE (blacs_matrix%p(iprow,ipcol)%block) + + END DO + END DO + +#else + + blacs_block => blacs_matrix%p(source,source)%block + + DO iblock_row=1,nblock_row + + block_node => first_block_node(matrix=sparse_matrix,& + block_row=iblock_row) + + DO WHILE (ASSOCIATED(block_node)) + + CALL get_block_node(block_node=block_node,& + block_col=iblock_col,& + block=sparse_block) + + icol = 1 + + DO icol_global=first_col(iblock_col),last_col(iblock_col) + + irow = 1 + + DO irow_global=first_row(iblock_row),last_row(iblock_row) + + sparse_block(irow,icol) = blacs_block(irow_global,icol_global) + + irow = irow + 1 + + END DO + + icol = icol + 1 + + END DO + + block_node => next_block_node(block_node) + + END DO + + END DO + +#endif + CALL timestop(0.0_wp,handle) + + END SUBROUTINE copy_blacs_to_sparse_matrix + +! ***************************************************************************** + + SUBROUTINE copy_sparse_to_blacs_matrix(sparse_matrix,blacs_matrix,context,& + globenv) + +! Purpose: Copy a sparse matrix to a BLACS matrix. The BLACS matrix blocks +! are allocated during the copy procedure. + +! History: - Creation (05.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(OUT) :: blacs_matrix + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(real_matrix_type), POINTER :: sparse_matrix + INTEGER, INTENT(IN) :: context + +! *** Local variables *** + + TYPE(real_block_node_type), POINTER :: block_node + + INTEGER :: group,handle,iblock_col,iblock_row,icol,icol_global,icol_local,& + ipcol,ipe,iprow,irow,irow_global,irow_local,jpcol,jprow,mypcol,& + mype,myprow,nblock_row,ncol_block,npcol,npe,nprow,nrow_block,& + source + + INTEGER, DIMENSION(:), POINTER :: first_col,first_row,last_col,last_row + REAL(wp), DIMENSION(:,:), POINTER :: blacs_block,sparse_block + +#if defined(__parallel) + INTEGER, EXTERNAL :: blacs_pnum,indxg2l,indxg2p + +#endif +! --------------------------------------------------------------------------- + + CALL timeset("copy_sparse_to_blacs_matrix","I","",handle) + + group = globenv%group + source = globenv%source + + CALL get_matrix_info(matrix=sparse_matrix,& + nblock_row=nblock_row,& + first_row=first_row,& + first_col=first_col,& + last_row=last_row,& + last_col=last_col) +#if defined(__parallel) + + CALL blacs_pinfo(mype,npe) + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + + nrow_block = blacs_matrix%nrow_block + ncol_block = blacs_matrix%ncol_block + + DO iprow=0,nprow-1 + DO ipcol=0,npcol-1 + + ipe = blacs_pnum(context,iprow,ipcol) + + IF (ipe /= mype) THEN + blacs_matrix%p(iprow,ipcol)%block =>& + reallocate(blacs_matrix%p(iprow,ipcol)%block,& + 1,blacs_matrix%p(iprow,ipcol)%nrow_local,& + 1,blacs_matrix%p(iprow,ipcol)%ncol_local) + ELSE + blacs_matrix%p(iprow,ipcol)%block(:,:) = 0.0_wp + END IF + + blacs_block => blacs_matrix%p(iprow,ipcol)%block + + DO iblock_row=1,nblock_row + + block_node => first_block_node(matrix=sparse_matrix,& + block_row=iblock_row) + + DO WHILE (ASSOCIATED(block_node)) + + CALL get_block_node(block_node=block_node,& + block_col=iblock_col,& + block=sparse_block) + + icol = 1 + + DO icol_global=first_col(iblock_col),last_col(iblock_col) + + jpcol = indxg2p(icol_global,ncol_block,mypcol,source,npcol) + + IF (jpcol == ipcol) THEN + + icol_local = indxg2l(icol_global,ncol_block,mypcol,source,& + npcol) + + irow = 1 + + DO irow_global=first_row(iblock_row),last_row(iblock_row) + + jprow = indxg2p(irow_global,nrow_block,myprow,source,nprow) + + IF (jprow == iprow) THEN + + irow_local = indxg2l(irow_global,nrow_block,myprow,source,& + nprow) + + blacs_block(irow_local,icol_local) = sparse_block(irow,& + icol) + + END IF + + irow = irow + 1 + + END DO + + END IF + + icol = icol + 1 + + END DO + + block_node => next_block_node(block_node) + + END DO + + END DO + + CALL mp_sum(blacs_block,ipe,group) + + IF (ipe /= mype) DEALLOCATE (blacs_matrix%p(iprow,ipcol)%block) + + END DO + END DO + +#else + + IF (.NOT.ASSOCIATED(blacs_matrix%p(source,source)%block)) THEN + blacs_matrix%p(source,source)%block =>& + reallocate(blacs_matrix%p(source,source)%block,& + 1,blacs_matrix%p(source,source)%nrow_local,& + 1,blacs_matrix%p(source,source)%ncol_local) + END IF + + blacs_block => blacs_matrix%p(source,source)%block + + blacs_block(:,:) = 0.0_wp + + DO iblock_row=1,nblock_row + + block_node => first_block_node(matrix=sparse_matrix,& + block_row=iblock_row) + + DO WHILE (ASSOCIATED(block_node)) + + CALL get_block_node(block_node=block_node,& + block_col=iblock_col,& + block=sparse_block) + + icol = 1 + + DO icol_global=first_col(iblock_col),last_col(iblock_col) + + irow = 1 + + DO irow_global=first_row(iblock_row),last_row(iblock_row) + + blacs_block(irow_global,icol_global) = sparse_block(irow,icol) + + irow = irow + 1 + + END DO + + icol = icol + 1 + + END DO + + block_node => next_block_node(block_node) + + END DO + + END DO + +#endif + CALL timestop(0.0_wp,handle) + + END SUBROUTINE copy_sparse_to_blacs_matrix + +! ***************************************************************************** + + SUBROUTINE deallocate_blacs_matrix(matrix) + +! Purpose: Deallocate a distributed BLACS matrix. + +! History: - Creation (08.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(INOUT) :: matrix + +! *** Local variables *** + + INTEGER :: ipcol,iprow,npcol,nprow + +! --------------------------------------------------------------------------- + + matrix%name = "" + + matrix%context = 0 + matrix%nrow_block = 0 + matrix%ncol_block = 0 + matrix%nrow_global = 0 + matrix%ncol_global = 0 + + IF (ASSOCIATED(matrix%descriptor)) DEALLOCATE (matrix%descriptor) + + IF (ASSOCIATED(matrix%p)) THEN + + nprow = SIZE(matrix%p,1) + npcol = SIZE(matrix%p,2) + + DO iprow=0,nprow-1 + DO ipcol=0,npcol-1 + IF (ASSOCIATED(matrix%p(iprow,ipcol)%block)) THEN + DEALLOCATE (matrix%p(iprow,ipcol)%block) + END IF + END DO + END DO + + DEALLOCATE (matrix%p) + + END IF + + END SUBROUTINE deallocate_blacs_matrix + +! ***************************************************************************** + + SUBROUTINE finish_blacs(context,globenv) + +! Purpose: Release the resources of a BLACS context. + +! History: - Creation (22.05.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(global_environment_type), INTENT(IN) :: globenv + INTEGER, INTENT(IN) :: context + +! *** Local variables *** + + INTEGER :: group,ipe,mype,npe,output_unit + LOGICAL :: ionode + + INTEGER, DIMENSION(:), POINTER :: pcon + +! --------------------------------------------------------------------------- + + group = globenv%group + ionode = globenv%ionode + output_unit = globenv%scr +#if defined(__parallel) + + IF (globenv%print%blacs_info) THEN + + CALL blacs_pinfo(mype,npe) + + NULLIFY (pcon) + pcon => reallocate(pcon,0,npe-1) + + pcon(mype) = context + + CALL mp_sum(pcon,group) + + IF (ionode) THEN + WRITE (UNIT=output_unit,FMT="(/,/,T2,A)")& + "BLACS INFORMATION (BLACS finished)" + WRITE (UNIT=output_unit,FMT="(/,T3,A,/)")& + " PE BLACS context" + WRITE (UNIT=output_unit,FMT="(I5,T10,I12)")& + (ipe,pcon(ipe),ipe=0,npe-1) + END IF + END IF + + CALL blacs_gridexit(context) + +#endif + END SUBROUTINE finish_blacs + +! ***************************************************************************** + + SUBROUTINE get_blacs_info(context,globenv,my_process_row,my_process_column,& + my_process_number,number_of_process_rows,& + number_of_process_columns,number_of_processes) + +! Purpose: Return informations about the specified BLACS context. + +! History: - Creation (19.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(global_environment_type), INTENT(IN) :: globenv + INTEGER, INTENT(IN) :: context + INTEGER, OPTIONAL, INTENT(OUT) :: my_process_column,& + my_process_number,& + my_process_row,& + number_of_process_columns,& + number_of_process_rows,& + number_of_processes + +! *** Local variables *** + + INTEGER :: mypcol,mype,myprow,npcol,npe,nprow,source + +#if defined(__parallel) + INTEGER, EXTERNAL :: blacs_pnum + +#endif +! --------------------------------------------------------------------------- + + source = globenv%source +#if defined(__parallel) + + CALL blacs_pinfo(mype,npe) + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + +#else + + myprow = source + mypcol = source + mype = source + nprow = 1 + npcol = 1 + npe = 1 + +#endif + IF (PRESENT(my_process_row)) my_process_row = myprow + IF (PRESENT(my_process_column)) my_process_column = mypcol + IF (PRESENT(my_process_number)) my_process_number = mype + IF (PRESENT(number_of_process_rows)) number_of_process_rows = nprow + IF (PRESENT(number_of_process_columns)) number_of_process_columns = npcol + IF (PRESENT(number_of_processes)) number_of_processes = npe + + END SUBROUTINE get_blacs_info + +! ***************************************************************************** + + SUBROUTINE get_blacs_matrix_info(matrix,name,nrow_global,ncol_global,& + nrow_block,ncol_block) + +! Purpose: Return informations about the specified BLACS matrix. + +! History: - Creation (08.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(IN) :: matrix + CHARACTER(LEN=*), OPTIONAL, INTENT(OUT) :: name + INTEGER, OPTIONAL, INTENT(OUT) :: ncol_block,ncol_global,& + nrow_block,nrow_global + +! --------------------------------------------------------------------------- + + IF (PRESENT(name)) name = matrix%name + IF (PRESENT(nrow_global)) nrow_global = matrix%nrow_global + IF (PRESENT(ncol_global)) ncol_global = matrix%ncol_global + IF (PRESENT(nrow_block)) nrow_block = matrix%nrow_block + IF (PRESENT(ncol_block)) ncol_block = matrix%ncol_block + + END SUBROUTINE get_blacs_matrix_info + +! ***************************************************************************** + + SUBROUTINE power_blacs_matrix(matrix,work,exponent,threshold,n_dependent,& + work_syevx,context,globenv) + +! Purpose: Raise the real symmetric n by n matrix to the power given by +! exponent. All eigenvectors with a corresponding eigenvalue lower +! than threshold are quenched. + +! History: - Creation (29.03.1999, Matthias Krack) +! - Parallelised using BLACS and ScaLAPACK (06.06.2001, MK) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(INOUT) :: matrix,work + TYPE(global_environment_type), INTENT(IN) :: globenv + REAL(wp), INTENT(IN) :: exponent,threshold,work_syevx + INTEGER, INTENT(IN) :: context + INTEGER, INTENT(OUT) :: n_dependent + +! *** Local variables *** + + REAL(wp) :: f,p + INTEGER :: handle,icol_global,icol_local,ipcol,iprow,irow_global,& + irow_local,mypcol,myprow,ncol_block,ncol_global,npcol,& + nprow,nrow_block,nrow_global,source + + REAL(wp), DIMENSION(:), POINTER :: eigenvalues + REAL(wp), DIMENSION(:,:), POINTER :: eigenvectors + +#if defined(__parallel) + INTEGER, EXTERNAL :: indxg2l,indxg2p + +#endif +! --------------------------------------------------------------------------- + + CALL timeset("power_blacs_matrix","I","",handle) + + source = globenv%source + n_dependent = 0 + p = 0.5_wp*exponent + + nrow_global = matrix%nrow_global + ncol_global = matrix%ncol_global + + NULLIFY (eigenvalues) + eigenvalues => reallocate(eigenvalues,1,ncol_global) + +! *** Compute the eigenvectors and eigenvalues *** + + CALL blacs_syevx(matrix,work,eigenvalues,ncol_global,work_syevx,context,& + globenv) +#if defined(__parallel) + + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + + nrow_block = work%nrow_block + ncol_block = work%ncol_block + + eigenvectors => work%p(myprow,mypcol)%block + +! *** Build matrix**exponent with eigenvector quenching *** + + p = 0.5_wp*exponent + + n_dependent = 0 + + DO icol_global=1,ncol_global + + IF (eigenvalues(icol_global) < threshold) THEN + + n_dependent = n_dependent + 1 + + ipcol = indxg2p(icol_global,ncol_block,mypcol,source,npcol) + + IF (mypcol == ipcol) THEN + icol_local = indxg2l(icol_global,ncol_block,mypcol,source,npcol) + DO irow_global=1,nrow_global + iprow = indxg2p(irow_global,nrow_block,myprow,source,nprow) + IF (myprow == iprow) THEN + irow_local = indxg2l(irow_global,nrow_block,myprow,source,nprow) + eigenvectors(irow_local,icol_local) = 0.0_wp + END IF + END DO + END IF + + ELSE + + f = eigenvalues(icol_global)**p + + ipcol = indxg2p(icol_global,ncol_block,mypcol,source,npcol) + + IF (mypcol == ipcol) THEN + icol_local = indxg2l(icol_global,ncol_block,mypcol,source,npcol) + DO irow_global=1,nrow_global + iprow = indxg2p(irow_global,nrow_block,myprow,source,nprow) + IF (myprow == iprow) THEN + irow_local = indxg2l(irow_global,nrow_block,myprow,source,nprow) + eigenvectors(irow_local,icol_local) =& + f*eigenvectors(irow_local,icol_local) + END IF + END DO + END IF + + END IF + + END DO + +#else + + eigenvectors => work%p(source,source)%block + +! *** Build matrix**exponent with eigenvector quenching *** + + DO icol_global=1,ncol_global + + IF (eigenvalues(icol_global) < threshold) THEN + + n_dependent = n_dependent + 1 + eigenvectors(1:nrow_global,icol_global) = 0.0_wp + + ELSE + + f = eigenvalues(icol_global)**p + eigenvectors(1:nrow_global,icol_global) =& + f*eigenvectors(1:nrow_global,icol_global) + + END IF + + END DO + +#endif + CALL blacs_syrk("U","N",ncol_global,1.0_wp,work,0.0_wp,matrix,context,& + globenv) + + DEALLOCATE (eigenvalues) + + CALL timestop(0.0_wp,handle) + + END SUBROUTINE power_blacs_matrix + +! ***************************************************************************** + + SUBROUTINE read_blacs_matrix(matrix,lunit,context,globenv) + +! Purpose: Read a BLACS matrix from the logical unit number "lunit". + +! History: - Creation (19.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(OUT) :: matrix + TYPE(global_environment_type), INTENT(IN) :: globenv + INTEGER, INTENT(IN) :: context,lunit + +! *** Local variables *** + + INTEGER :: i,j,mypcol,myprow,ncol_local,npcol,nprow,nrow_local,source + +! --------------------------------------------------------------------------- + + source = globenv%source +#if defined(__parallel) + + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + +#else + + myprow = source + mypcol = source + +#endif + nrow_local = matrix%p(myprow,mypcol)%nrow_local + ncol_local = matrix%p(myprow,mypcol)%ncol_local + + READ (UNIT=lunit) ((matrix%p(myprow,mypcol)%block(i,j),i=1,nrow_local),& + j=1,ncol_local) + + END SUBROUTINE read_blacs_matrix + +! ***************************************************************************** + + SUBROUTINE replicate_blacs_matrix(prototype_matrix,new_matrix,name) + +! Purpose: Allocate a distributed BLACS matrix using a prototype matrix. + +! History: - Creation (08.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(IN) :: prototype_matrix + TYPE(blacs_matrix_type), INTENT(OUT) :: new_matrix + CHARACTER(LEN=*), INTENT(IN) :: name + +! *** Local variables *** + + INTEGER :: ipcol,iprow,npcol,nprow + +! --------------------------------------------------------------------------- + + new_matrix%name = name + + new_matrix%context = prototype_matrix%context + new_matrix%nrow_block = prototype_matrix%nrow_block + new_matrix%ncol_block = prototype_matrix%ncol_block + new_matrix%nrow_global = prototype_matrix%nrow_global + new_matrix%ncol_global = prototype_matrix%ncol_global + NULLIFY (new_matrix%descriptor) + new_matrix%descriptor => reallocate(new_matrix%descriptor,1,9) + new_matrix%descriptor(:) = prototype_matrix%descriptor(:) + + nprow = SIZE(prototype_matrix%p,1) + npcol = SIZE(prototype_matrix%p,2) + + ALLOCATE (new_matrix%p(0:nprow-1,0:npcol-1)) + + DO iprow=0,nprow-1 + DO ipcol=0,npcol-1 + NULLIFY (new_matrix%p(iprow,ipcol)%block) + new_matrix%p(iprow,ipcol)%nrow_local =& + prototype_matrix%p(iprow,ipcol)%nrow_local + new_matrix%p(iprow,ipcol)%ncol_local =& + prototype_matrix%p(iprow,ipcol)%ncol_local + IF (ASSOCIATED(prototype_matrix%p(iprow,ipcol)%block)) THEN + new_matrix%p(iprow,ipcol)%block =>& + reallocate(new_matrix%p(iprow,ipcol)%block,& + 1,new_matrix%p(iprow,ipcol)%nrow_local,& + 1,new_matrix%p(iprow,ipcol)%ncol_local) + new_matrix%p(iprow,ipcol)%block(:,:) =& + prototype_matrix%p(iprow,ipcol)%block(:,:) + END IF + END DO + END DO + + END SUBROUTINE replicate_blacs_matrix + +! ***************************************************************************** + + SUBROUTINE start_blacs(nprow,npcol,context,globenv) + +! Purpose: Initialize a BLACS process grid. The BLACS context is returned. + +! History: - Creation (22.05.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(global_environment_type), INTENT(IN) :: globenv + INTEGER, INTENT(INOUT) :: context,npcol,nprow + +! *** Local variables *** + + INTEGER :: group,ipe,mypcol,mype,myprow,npe,output_unit,source + LOGICAL :: ionode + + INTEGER, DIMENSION(:), POINTER :: pcol,pcon,prow + +! --------------------------------------------------------------------------- + + group = globenv%group + ionode = globenv%ionode + output_unit = globenv%scr + source = globenv%source +#if defined(__parallel) + + CALL blacs_pinfo(mype,npe) + CALL blacs_get(-1,0,context) + + IF (nprow*npcol /= npe) THEN + DO ipe=CEILING(SQRT(REAL(npe,wp))),npe + IF (MODULO(npe,ipe) == 0) THEN + nprow = ipe + npcol = npe/nprow + EXIT + END IF + END DO + END IF + + CALL blacs_gridinit(context,"Row-major",nprow,npcol) + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + + IF (globenv%print%blacs_info) THEN + + NULLIFY (prow,pcol,pcon) + prow => reallocate(prow,0,npe-1) + pcol => reallocate(pcol,0,npe-1) + pcon => reallocate(pcon,0,npe-1) + + prow(mype) = myprow + pcol(mype) = mypcol + pcon(mype) = context + + CALL mp_sum(prow,group) + CALL mp_sum(pcol,group) + CALL mp_sum(pcon,group) + + IF (ionode) THEN + WRITE (UNIT=output_unit,FMT="(/,/,T2,A,/,/,(T3,A,T32,I6))")& + "BLACS INFORMATION (BLACS started)",& + "Number of processes: ",nprow*npcol,& + "Number of process rows: ",nprow,& + "Number of process columns: ",npcol + WRITE (UNIT=output_unit,FMT="(/,T3,A,/)")& + " PE process row process column BLACS context" + WRITE (UNIT=output_unit,FMT="(I5,T14,I6,T31,I6,T41,I12)")& + (ipe,prow(ipe),pcol(ipe),pcon(ipe),ipe=0,npe-1) + END IF + + DEALLOCATE (prow,pcol,pcon) + + END IF + +#endif + END SUBROUTINE start_blacs + +! ***************************************************************************** + + SUBROUTINE symmetrise_blacs_matrix(matrix,work,context,globenv) + +! Purpose: Symmetrise a symmetric BLACS matrix. + +! History: - Creation (12.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(INOUT) :: matrix,work + TYPE(global_environment_type), INTENT(IN) :: globenv + INTEGER, INTENT(IN) :: context + +! *** Local variables *** + + INTEGER :: handle,icol_global,icol_local,ipcol,iprow,irow_global,& + irow_local,mypcol,myprow,ncol_block,ncol_global,ncol_local,& + npcol,nprow,nrow_block,nrow_global,nrow_local,source + + INTEGER, DIMENSION(:), POINTER :: desca,descc + REAL(wp), DIMENSION(:,:), POINTER :: a,c + +#if defined(__parallel) + INTEGER, EXTERNAL :: indxl2g + +#endif +! --------------------------------------------------------------------------- + + CALL timeset("symmetrise_blacs_matrix","I","",handle) + + source = globenv%source +#if defined(__parallel) + + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + + nrow_global = matrix%nrow_global + ncol_global = matrix%ncol_global + + nrow_block = matrix%nrow_block + ncol_block = matrix%ncol_block + + nrow_local = matrix%p(myprow,mypcol)%nrow_local + ncol_local = matrix%p(myprow,mypcol)%ncol_local + + a => work%p(myprow,mypcol)%block + desca => work%descriptor + c => matrix%p(myprow,mypcol)%block + descc => matrix%descriptor + + DO icol_local=1,ncol_local + icol_global = indxl2g(icol_local,ncol_block,mypcol,source,npcol) + DO irow_local=1,nrow_local + irow_global = indxl2g(irow_local,nrow_block,myprow,source,nprow) + IF (irow_global > icol_global) THEN + c(irow_local,icol_local) = 0.0_wp + ELSE IF (irow_global == icol_global) THEN + c(irow_local,icol_local) = 0.5_wp*c(irow_local,icol_local) + END IF + END DO + END DO + + a(:,:) = c(:,:) + + CALL pdtran(nrow_global,ncol_global,1.0_wp,a,1,1,desca,1.0_wp,c,1,1,descc) + +#else + + a => matrix%p(source,source)%block + + CALL symmetrize_matrix(a,"upper_to_lower") + +#endif + CALL timestop(0.0_wp,handle) + + END SUBROUTINE symmetrise_blacs_matrix + +! ***************************************************************************** + + SUBROUTINE write_blacs_matrix(matrix,lunit,context,globenv) + +! Purpose: Write a BLACS matrix to the logical unit number "lunit". + +! History: - Creation (19.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(IN) :: matrix + TYPE(global_environment_type), INTENT(IN) :: globenv + INTEGER, INTENT(IN) :: context,lunit + +! *** Local variables *** + + INTEGER :: i,j,mypcol,myprow,ncol_local,npcol,nprow,nrow_local,source + +! --------------------------------------------------------------------------- + + source = globenv%source +#if defined(__parallel) + + CALL blacs_gridinfo(context,nprow,npcol,myprow,mypcol) + +#else + + myprow = source + mypcol = source + +#endif + nrow_local = matrix%p(myprow,mypcol)%nrow_local + ncol_local = matrix%p(myprow,mypcol)%ncol_local + + WRITE (UNIT=lunit) ((matrix%p(myprow,mypcol)%block(i,j),i=1,nrow_local),& + j=1,ncol_local) + + END SUBROUTINE write_blacs_matrix + +! ***************************************************************************** + +END MODULE blacs diff --git a/src/collocate_density.F b/src/collocate_density.F index 77cc073422..a9c30986a6 100644 --- a/src/collocate_density.F +++ b/src/collocate_density.F @@ -30,10 +30,12 @@ MODULE collocate_density USE atoms, ONLY: atom_info USE basis_set_types, ONLY: gto_basis_set_type,maxco,maxsgf,maxsgf_set USE coefficient_types, ONLY: coeff_type + USE global_types, ONLY: global_environment_type USE interactions, ONLY: eps_rho_gspace,eps_rho_rspace,exp_radius USE matrix_types, ONLY: get_block_node,& real_matrix_set_type USE memory_utilities, ONLY: reallocate + USE message_passing, ONLY: mp_sum USE neighbor_list_types, ONLY: extract_neighbor_list,& find_neighbor_list,& first_neighbor_list,& @@ -70,10 +72,11 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE calculate_rho_core(rho,total_rho) + SUBROUTINE calculate_rho_core(rho,total_rho,globenv) - TYPE(coeff_type), INTENT(INOUT) :: rho - REAL(wp), INTENT(OUT) :: total_rho + TYPE(coeff_type), INTENT(INOUT) :: rho + TYPE(global_environment_type), INTENT(IN) :: globenv + REAL(wp), INTENT(OUT) :: total_rho ! *** Local variables *** @@ -110,6 +113,12 @@ CONTAINS END DO +! IF (ASSOCIATED(rho%pw%cc3d)) THEN +! CALL mp_sum(rho%pw%cc3d,globenv%group) +! ELSE IF (ASSOCIATED(rho%pw%cr3d)) THEN +! CALL mp_sum(rho%pw%cr3d,globenv%group) +! END IF + total_rho = calculate_total_rho(rho) CALL timestop(0.0_wp,handle) @@ -118,12 +127,13 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE calculate_rho_elec(p,rho,total_rho,nproduct) + SUBROUTINE calculate_rho_elec(p,rho,total_rho,nproduct,globenv) - TYPE(real_matrix_set_type), INTENT(IN) :: p - TYPE(coeff_type), INTENT(INOUT) :: rho - REAL(wp), INTENT(OUT) :: total_rho - INTEGER, OPTIONAL, INTENT(OUT) :: nproduct + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(coeff_type), INTENT(INOUT) :: rho + TYPE(real_matrix_set_type), INTENT(IN) :: p + REAL(wp), INTENT(OUT) :: total_rho + INTEGER, INTENT(OUT) :: nproduct ! *** Local variables *** @@ -143,7 +153,7 @@ CONTAINS INTEGER, DIMENSION(:), POINTER :: lb_grid,orb_neighbors REAL(wp), DIMENSION(:), POINTER :: orb_r2,rpgfa,rpgfb,zeta,zetb - REAL(wp), DIMENSION(:,:), POINTER :: orb_r,p_block,pab,work + REAL(wp), DIMENSION(:,:), POINTER :: orb_r,p_block,pab,s_block,work ! --------------------------------------------------------------------------- @@ -201,6 +211,8 @@ CONTAINS block_col=jatom,& block=p_block) + IF (.NOT.ASSOCIATED(p_block)) CYCLE + IF (MAXVAL(ABS(p_block(:,:))) < 1.0E-14_wp) CYCLE DO iset=1,nseta @@ -311,12 +323,22 @@ CONTAINS END DO - total_rho = calculate_total_rho(rho) - - IF (PRESENT(nproduct)) nproduct = npgf_product +! *** Release work storage *** DEALLOCATE (pab,work) + CALL mp_sum(npgf_product,globenv%group) + + nproduct = npgf_product + + IF (ASSOCIATED(rho%pw%cc3d)) THEN + CALL mp_sum(rho%pw%cc3d,globenv%group) + ELSE IF (ASSOCIATED(rho%pw%cr3d)) THEN + CALL mp_sum(rho%pw%cr3d,globenv%group) + END IF + + total_rho = calculate_total_rho(rho) + CALL timestop(0.0_wp,handle) END SUBROUTINE calculate_rho_elec @@ -325,7 +347,7 @@ CONTAINS FUNCTION calculate_total_rho(rho) RESULT(total_rho) - TYPE(coeff_type), INTENT(IN), TARGET :: rho + TYPE(coeff_type), TARGET, INTENT(IN) :: rho REAL(wp) :: total_rho @@ -357,9 +379,9 @@ CONTAINS USE mathconstants, ONLY: pi,twopi TYPE(pw_type), TARGET, INTENT(INOUT) :: pw - REAL(wp), INTENT(IN) :: rab2,scale,zeta,zetb - INTEGER, INTENT(IN) :: la_max,la_min,lb_max,lb_min - REAL(wp), DIMENSION(3), INTENT(IN) :: ra,rab + REAL(wp), INTENT(IN) :: rab2,scale,zeta,zetb + INTEGER, INTENT(IN) :: la_max,la_min,lb_max,lb_min + REAL(wp), DIMENSION(3), INTENT(IN) :: ra,rab REAL(wp), DIMENSION(ncoset(la_max),ncoset(lb_max)), INTENT(IN) :: pab @@ -544,9 +566,9 @@ CONTAINS ra,rab,rab2,scale,pab,pw) TYPE(pw_type), TARGET, INTENT(INOUT) :: pw - REAL(wp), INTENT(IN) :: rab2,scale,zeta,zetb - INTEGER, INTENT(IN) :: la_max,la_min,lb_max,lb_min - REAL(wp), DIMENSION(3), INTENT(IN) :: ra,rab + REAL(wp), INTENT(IN) :: rab2,scale,zeta,zetb + INTEGER, INTENT(IN) :: la_max,la_min,lb_max,lb_min + REAL(wp), DIMENSION(3), INTENT(IN) :: ra,rab REAL(wp), DIMENSION(ncoset(la_max),ncoset(lb_max)), INTENT(IN) :: pab diff --git a/src/core_energies.F b/src/core_energies.F index 7f92db5033..62f85a4063 100644 --- a/src/core_energies.F +++ b/src/core_energies.F @@ -25,6 +25,9 @@ MODULE core_energies USE kinds, ONLY: wp => dp + USE global_types, ONLY: global_environment_type + USE message_passing, ONLY: mp_sum + IMPLICIT NONE PRIVATE @@ -39,11 +42,13 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE calculate_ecore(h,p,ecore) + SUBROUTINE calculate_ecore(h,p,ecore,globenv) ! Purpose: Calculate the core Hamiltonian energy which includes the kinetic ! and the potential energy of the electrons. +! History: - Creation (03.05.2001, Matthias Krack) + ! *************************************************************************** USE matrix_types, ONLY: first_block_node,& @@ -53,8 +58,9 @@ CONTAINS real_block_node_type,& real_matrix_set_type - TYPE(real_matrix_set_type), INTENT(IN) :: h,p - REAL(wp), INTENT(OUT) :: ecore + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(real_matrix_set_type), INTENT(IN) :: h,p + REAL(wp), INTENT(OUT) :: ecore ! *** Local variables *** @@ -105,21 +111,26 @@ CONTAINS END DO + CALL mp_sum(ecore,globenv%group) + END SUBROUTINE calculate_ecore ! ***************************************************************************** - SUBROUTINE calculate_ecore_overlap(ecore_overlap) + SUBROUTINE calculate_ecore_overlap(ecore_overlap,globenv) ! Purpose: Calculate the overlap energy of the core charge distribution. +! History: - Creation (30.04.2001, Matthias Krack) + ! *************************************************************************** USE atomic_kinds, ONLY: kind_info USE atoms, ONLY: atom_info,natom USE cell_parameters, ONLY: abc - REAL(wp), INTENT(OUT) :: ecore_overlap + TYPE(global_environment_type), INTENT(IN) :: globenv + REAL(wp), INTENT(OUT) :: ecore_overlap ! *** Local variables *** @@ -137,6 +148,8 @@ CONTAINS DO iatom=1,natom + IF (globenv%mepos /= MODULO(iatom,globenv%num_pe)) CYCLE + ikind = atom_info(iatom)%kind alpha_a = kind_info(ikind)%alpha_core_charge @@ -186,6 +199,8 @@ CONTAINS END DO + CALL mp_sum(ecore_overlap,globenv%group) + END SUBROUTINE calculate_ecore_overlap ! ***************************************************************************** @@ -194,6 +209,8 @@ CONTAINS ! Purpose: Calculate the self energy of the core charge distribution. +! History: - Creation (27.04.2001, Matthias Krack) + ! *************************************************************************** USE atomic_kinds, ONLY: kind_info,nkind diff --git a/src/core_hamiltonian.F b/src/core_hamiltonian.F index d72d44ea5e..4953f1ded2 100644 --- a/src/core_hamiltonian.F +++ b/src/core_hamiltonian.F @@ -102,11 +102,13 @@ MODULE core_hamiltonian deallocate_matrix_row,& first_block_node,& get_block_node,& + get_matrix_info,& next_block_node,& real_block_node_type,& real_matrix_set_type,& replicate_matrix_structure USE memory_utilities, ONLY: reallocate + USE message_passing, ONLY: mp_sum USE method_specifications, ONLY: allchem,gpw,maxder USE neighbor_list_types, ONLY: extract_neighbor_list,& find_neighbor_list,& @@ -286,14 +288,23 @@ CONTAINS TYPE(neighbor_list_type), POINTER :: neighbor_list TYPE(neighbor_node_type), POINTER :: neighbor_node - INTEGER :: iatom,ij,ikind,ipe,istat,jatom,jkind,jpe,nsgfa,nsgfb + INTEGER :: group,iatom,ij,ikind,ipe,istat,jatom,jkind,jpe,mype,& + nblock_full_matrix,npe,nsgfa,nsgfb,output_unit,sum_nblock_pe + LOGICAL :: ionode - INTEGER, DIMENSION(natom*(natom+1)/2) :: nelement_block - INTEGER, DIMENSION(0:globenv%num_pe-1) :: nelement_pe + INTEGER, DIMENSION(natom*(natom+1)/2) :: nblock,nelement + INTEGER, DIMENSION(0:globenv%num_pe-1) :: nblock_pe,nelement_pe ! --------------------------------------------------------------------------- - nelement_block(:) = 0 + group = globenv%group + ionode = globenv%ionode + mype = globenv%mepos + npe = globenv%num_pe + output_unit = globenv%scr + + nblock(:) = 0 + nelement(:) = 0 ! *** Allocate the overlap matrix *** @@ -335,7 +346,8 @@ CONTAINS IF (iatom <= jatom) THEN ij = iatom + jatom*(jatom - 1)/2 - nelement_block(ij) = nelement_block(ij) + nsgfa*nsgfb + nblock(ij) = nblock(ij) + 1 + nelement(ij) = nelement(ij) + nsgfa*nsgfb END IF neighbor_node => next_neighbor_node(neighbor_node) @@ -348,6 +360,7 @@ CONTAINS ! *** Distribute the atom blocks *** + nblock_pe(:) = 0 nelement_pe(:) = 0 DO iatom=1,natom @@ -355,31 +368,68 @@ CONTAINS ij = iatom + jatom*(jatom - 1)/2 - IF (nelement_block(ij) > 0) THEN + IF (nelement(ij) > 0) THEN ipe = 0 - DO jpe=0,globenv%num_pe-1 + DO jpe=0,npe-1 IF (nelement_pe(jpe) < nelement_pe(ipe)) ipe = jpe END DO - nelement_pe(ipe) = nelement_pe(ipe) + nelement_block(ij) + nblock_pe(ipe) = nblock_pe(ipe) + nblock(ij) + nelement_pe(ipe) = nelement_pe(ipe) + nelement(ij) - IF (ipe == globenv%mepos) CALL add_block_node(matrix=s%matrix,& - block_row=iatom,& - block_col=jatom) + IF (ipe == mype) CALL add_block_node(matrix=s%matrix,& + block_row=iatom,& + block_col=jatom) END IF END DO END DO -! *** Print the distribution *** +! *** Print the distribution of the overlap matrix *** + + IF (globenv%print%distribution) THEN + + IF (ionode) THEN + WRITE (UNIT=output_unit,& + FMT="(/,/,T2,A,/,/,T3,A,/,/,T5,A,/,/,(I6,8X,I8,8X,I10))")& + "DISTRIBUTION OF THE OVERLAP MATRIX ELEMENTS",& + "Image atoms included:",& + "PE Matrix blocks Matrix elements",& + (ipe,nblock_pe(ipe),nelement_pe(ipe),ipe=0,npe-1) + WRITE (UNIT=output_unit,FMT="(/,T4,A3,8X,I8,8X,I10)")& + "Sum",SUM(nblock_pe),SUM(nelement_pe) + END IF + + nblock_pe(:) = 0 + nelement_pe(:) = 0 + + CALL get_matrix_info(matrix=s%matrix,& + nblock_allocated=nblock_pe(mype),& + nelement_allocated=nelement_pe(mype)) + + CALL mp_sum(nblock_pe,group) + CALL mp_sum(nelement_pe,group) + + nblock_full_matrix = natom*(natom + 1)/2 + sum_nblock_pe = SUM(nblock_pe) + + IF (ionode) THEN + WRITE (UNIT=output_unit,& + FMT="(/,/,T3,A,/,/,T5,A,/,/,(I6,8X,I8,8X,I10))")& + "Image atoms not included:",& + "PE Matrix blocks Matrix elements",& + (ipe,nblock_pe(ipe),nelement_pe(ipe),ipe=0,npe-1) + WRITE (UNIT=output_unit,FMT="(/,T4,A3,8X,I8,8X,I10)")& + "Sum",sum_nblock_pe,SUM(nelement_pe) + WRITE (UNIT=output_unit,FMT="(/,T4,A3,8X,I8,A,F5.1,A)")& + " of",nblock_full_matrix," blocks in the full matrix (",& + 100.0_wp*REAL(sum_nblock_pe,wp)/REAL(nblock_full_matrix,wp),& + " % occupation)" + END IF - IF (globenv%ionode.AND.globenv%print%distribution) THEN - WRITE (UNIT=globenv%scr,FMT="(/,T5,A,/,/,(I6,8X,I10))")& - "PE Overlap matrix elements",& - (ipe,nelement_pe(ipe),ipe=0,globenv%num_pe-1) END IF END SUBROUTINE distribute_s_matrix diff --git a/src/diis.F b/src/diis.F index dd0f3977b2..43f40eee4c 100644 --- a/src/diis.F +++ b/src/diis.F @@ -14,7 +14,7 @@ !! Matthias Krack (28.06.2000) !! !! MODIFICATION HISTORY -!! none +!! Changed to BLACS matrix usage (08.06.2001, MK) !! !! SOURCE !****************************************************************************** @@ -23,18 +23,24 @@ MODULE diis USE kinds, ONLY: wp => dp - USE core_hamiltonian, ONLY: s - USE global_types, ONLY: global_environment_type - USE mathlib, ONLY: diagonalize_matrix,& - symmetrize_matrix - USE matrix_types, ONLY: allocate_matrix,& - copy_matrix,& - get_block_node,& - get_matrix_info,& - real_matrix_set_type - USE memory_utilities, ONLY: reallocate - USE mo_types, ONLY: mo_set_type - USE timings, ONLY: timeset,timestop + USE blacs, ONLY: blacs_add,& + blacs_gemm,& + blacs_matrix_type,& + blacs_maxval,& + blacs_set_all,& + blacs_symm,& + blacs_trace,& + copy_blacs_to_blacs_matrix,& + copy_sparse_to_blacs_matrix,& + get_blacs_matrix_info,& + replicate_blacs_matrix + USE core_hamiltonian, ONLY: s + USE global_types, ONLY: global_environment_type + USE mathlib, ONLY: diagonalize_matrix,& + symmetrize_matrix + USE memory_utilities, ONLY: reallocate + USE mo_types, ONLY: mo_set_type + USE timings, ONLY: timeset,timestop IMPLICIT NONE @@ -42,11 +48,13 @@ MODULE diis TYPE diis_buffer_type PRIVATE - INTEGER :: nbuffer,ncall - TYPE(real_matrix_set_type), DIMENSION(:), POINTER :: error,parameter - REAL(wp), DIMENSION(:,:), POINTER :: b_matrix + INTEGER :: nbuffer,ncall + TYPE(blacs_matrix_type), DIMENSION(:), POINTER :: error,parameter + REAL(wp), DIMENSION(:,:), POINTER :: b_matrix END TYPE diis_buffer_type + TYPE(blacs_matrix_type), POINTER :: new_errors,old_errors,parameters + TYPE(diis_buffer_type) :: scf_diis_buffer REAL(wp) :: eps_diis = 0.1_wp INTEGER :: max_diis = 4 @@ -66,16 +74,18 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE allocate_scf_diis_buffer(nbuffer,nao) + SUBROUTINE allocate_scf_diis_buffer(nbuffer,nao,prototype_matrix) ! Purpose: Allocate and initialize a DIIS buffer for "nao*nao" parameter ! variables and with a buffer size of "nbuffer". ! History: - Creation (07.05.2001, Matthias Krack) +! - Changed to BLACS matrix usage (08.06.2001, MK) ! *************************************************************************** - INTEGER, INTENT(IN) :: nbuffer,nao + TYPE(blacs_matrix_type), INTENT(IN) :: prototype_matrix + INTEGER, INTENT(IN) :: nbuffer,nao ! *** Local variables *** @@ -90,18 +100,14 @@ CONTAINS ALLOCATE (scf_diis_buffer%parameter(nbuffer)) DO ibuffer=1,nbuffer - NULLIFY (scf_diis_buffer%error(ibuffer)%matrix) - CALL allocate_matrix(matrix=scf_diis_buffer%error(ibuffer)%matrix,& - nrow=nao,& - ncol=nao,& - matrix_name="SCF DIIS ERROR MATRIX",& - matrix_symmetry="symmetric") - NULLIFY (scf_diis_buffer%parameter(ibuffer)%matrix) - CALL allocate_matrix(matrix=scf_diis_buffer%parameter(ibuffer)%matrix,& - nrow=nao,& - ncol=nao,& - matrix_name="SCF DIIS PARAMETER MATRIX",& - matrix_symmetry="symmetric") + new_errors => scf_diis_buffer%error(ibuffer) + CALL replicate_blacs_matrix(prototype_matrix=prototype_matrix,& + new_matrix=new_errors,& + name="SCF DIIS ERROR MATRIX") + parameters => scf_diis_buffer%parameter(ibuffer) + CALL replicate_blacs_matrix(prototype_matrix=prototype_matrix,& + new_matrix=parameters,& + name="SCF DIIS PARAMETER MATRIX") END DO NULLIFY (scf_diis_buffer%b_matrix) @@ -115,29 +121,30 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE scf_diis(globenv,c,kc,sc,delta,error_max,diis_step) + SUBROUTINE scf_diis(c,kc,sc,delta,error_max,diis_step,context,globenv) ! Purpose: Update the SCF DIIS buffer. ! History: - Creation (07.05.2001, Matthias Krack) +! - Changed to BLACS matrix usage (08.06.2001, MK) ! *************************************************************************** TYPE(global_environment_type), INTENT(IN) :: globenv TYPE(mo_set_type), INTENT(IN) :: c - TYPE(real_matrix_set_type), INTENT(INOUT) :: kc,sc + TYPE(blacs_matrix_type), INTENT(INOUT) :: kc,sc REAL(wp), INTENT(IN) :: delta REAL(wp), INTENT(OUT) :: error_max + INTEGER, INTENT(IN) :: context LOGICAL, INTENT(OUT) :: diis_step ! *** Local variables *** INTEGER :: handle,ib,jb,nao,nb,nb1,nmo,nparameter,output_unit + LOGICAL :: ionode REAL(wp), DIMENSION(:), POINTER :: ev - REAL(wp), DIMENSION(:,:), POINTER :: a,b,c_matrix,kc_matrix,& - new_error_matrix,old_error_matrix,& - parameter_matrix,sc_matrix + REAL(wp), DIMENSION(:,:), POINTER :: a,b ! --------------------------------------------------------------------------- @@ -149,16 +156,18 @@ CONTAINS IF (max_diis < 1) RETURN + ionode = globenv%ionode output_unit= globenv%scr - CALL get_matrix_info(matrix=c%eigenvectors%matrix,nrow=nao) + CALL get_blacs_matrix_info(matrix=c%eigenvectors,& + nrow_global=nao) nmo = c%homo nparameter = nao*nao IF (.NOT.init_scf_diis_buffer_done) THEN - CALL allocate_scf_diis_buffer(max_diis,nao) - IF (globenv%print%diis_information) THEN + CALL allocate_scf_diis_buffer(max_diis,nao,kc) + IF (ionode.AND.globenv%print%diis_information) THEN WRITE (UNIT=output_unit,FMT="(/,T2,A)")& "The SCF DIIS buffer was allocated and initialised" END IF @@ -168,52 +177,30 @@ CONTAINS scf_diis_buffer%ncall = scf_diis_buffer%ncall + 1 nb = MIN(scf_diis_buffer%ncall,scf_diis_buffer%nbuffer) - CALL get_block_node(matrix=scf_diis_buffer%parameter(ib)%matrix,& - block_row=1,& - block_col=1,& - block=parameter_matrix) + parameters => scf_diis_buffer%parameter(ib) + CALL copy_blacs_to_blacs_matrix(kc,parameters) + CALL blacs_symm("L","U",nao,nmo,2.0_wp,parameters,c%eigenvectors,& + 0.0_wp,kc,context,globenv) +!MK CALL blacs_gemm("T","N",nao,nmo,nao,2.0_wp,parameters,c%eigenvectors,& +!MK 0.0_wp,kc,context,globenv) - CALL get_block_node(matrix=scf_diis_buffer%error(ib)%matrix,& - block_row=1,& - block_col=1,& - block=new_error_matrix) + new_errors => scf_diis_buffer%error(ib) + CALL copy_sparse_to_blacs_matrix(s%matrix,new_errors,context,globenv) + CALL blacs_symm("L","U",nao,nmo,2.0_wp,new_errors,c%eigenvectors,& + 0.0_wp,sc,context,globenv) +!MK CALL blacs_gemm("T","N",nao,nmo,nao,2.0_wp,new_errors,c%eigenvectors,& +!MK 0.0_wp,sc,context,globenv) - CALL get_block_node(matrix=c%eigenvectors%matrix,& - block_row=1,& - block_col=1,& - block=c_matrix) - - CALL get_block_node(matrix=kc%matrix,& - block_row=1,& - block_col=1,& - block=kc_matrix) - - CALL get_block_node(matrix=sc%matrix,& - block_row=1,& - block_col=1,& - block=sc_matrix) - - parameter_matrix(:,:) = kc_matrix(:,:) - - new_error_matrix(:,:) = 0.0_wp - CALL copy_matrix(s%matrix,scf_diis_buffer%error(ib)%matrix) - CALL symmetrize_matrix(new_error_matrix,"upper_to_lower") - - CALL dgemm("T","N",nao,nmo,nao,2.0_wp,parameter_matrix,nao,c_matrix,nao,& - 0.0_wp,kc_matrix,nao) - CALL dgemm("T","N",nao,nmo,nao,2.0_wp,new_error_matrix,nao,c_matrix,nao,& - 0.0_wp,sc_matrix,nao) - - CALL dgemm("N","T",nao,nao,nmo,1.0_wp,sc_matrix,nao,kc_matrix,nao,0.0_wp,& - new_error_matrix,nao) - CALL dgemm("N","T",nao,nao,nmo,1.0_wp,kc_matrix,nao,sc_matrix,nao,-1.0_wp,& - new_error_matrix,nao) + CALL blacs_gemm("N","T",nao,nao,nmo,1.0_wp,sc,kc, 0.0_wp,new_errors,& + context,globenv) + CALL blacs_gemm("N","T",nao,nao,nmo,1.0_wp,kc,sc,-1.0_wp,new_errors,& + context,globenv) ! *** Get maximum error *** - error_max = MAXVAL(ABS(new_error_matrix)) + CALL blacs_maxval(new_errors,error_max,context,globenv) - IF (globenv%print%diis_information) THEN + IF (ionode.AND.globenv%print%diis_information) THEN WRITE (UNIT=output_unit,FMT="(/,T2,A,E12.3)")& "Maximum SCF DIIS error vector element:",error_max END IF @@ -226,19 +213,15 @@ CONTAINS IF (error_max < eps_diis) THEN + b => scf_diis_buffer%b_matrix + DO jb=1,nb - CALL get_block_node(matrix=scf_diis_buffer%error(jb)%matrix,& - block_row=1,& - block_col=1,& - block=old_error_matrix) - CALL dgemm("T","N",1,1,nparameter,1.0_wp,old_error_matrix,nparameter,& - new_error_matrix,nparameter,0.0_wp,& - scf_diis_buffer%b_matrix(jb,ib),& - SIZE(scf_diis_buffer%b_matrix,1)) - scf_diis_buffer%b_matrix(ib,jb) = scf_diis_buffer%b_matrix(jb,ib) + old_errors => scf_diis_buffer%error(jb) + CALL blacs_trace(old_errors,new_errors,b(jb,ib),context,globenv) + b(ib,jb) = b(jb,ib) END DO - IF (globenv%print%diis_information) THEN + IF (ionode.AND.globenv%print%diis_information) THEN WRITE (UNIT=output_unit,FMT="(/,T2,A)")& "The SCF DIIS buffer was updated" END IF @@ -256,7 +239,6 @@ CONTAINS nb1 = nb + 1 NULLIFY (a,b,ev) - a => reallocate(a,1,nb1,1,nb1) b => reallocate(b,1,nb1,1,nb1) ev => reallocate(ev,1,nb1) @@ -287,21 +269,18 @@ CONTAINS ! *** Update Kohn-Sham matrix *** - kc_matrix(:,:) = 0.0_wp + CALL blacs_set_all(kc,0.0_wp,context,globenv) DO jb=1,nb - CALL get_block_node(matrix=scf_diis_buffer%parameter(jb)%matrix,& - block_row=1,& - block_col=1,& - block=parameter_matrix) - kc_matrix(:,:) = kc_matrix(:,:) - ev(jb)*parameter_matrix(:,:) + parameters => scf_diis_buffer%parameter(jb) + CALL blacs_add(1.0_wp,kc,-ev(jb),parameters,context,globenv) END DO DEALLOCATE (a,b,ev) ELSE - kc_matrix(:,:) = parameter_matrix(:,:) + CALL copy_blacs_to_blacs_matrix(parameters,kc) END IF diff --git a/src/integrate_potential.F b/src/integrate_potential.F index db336e0984..70f1e5d36d 100644 --- a/src/integrate_potential.F +++ b/src/integrate_potential.F @@ -30,10 +30,12 @@ MODULE integrate_potential USE atoms, ONLY: atom_info USE basis_set_types, ONLY: gto_basis_set_type,maxco,maxsgf,maxsgf_set USE coefficient_types, ONLY: coeff_type + USE global_types, ONLY: global_environment_type USE mathlib, ONLY: symmetrize_matrix USE matrix_types, ONLY: get_block_node,& real_matrix_set_type USE memory_utilities, ONLY: reallocate + USE message_passing, ONLY: mp_sum USE neighbor_list_types, ONLY: extract_neighbor_list,& find_neighbor_list,& first_neighbor_list,& @@ -63,11 +65,12 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE integrate_v_rspace(v_rspace,h,nproduct) + SUBROUTINE integrate_v_rspace(v_rspace,h,nproduct,globenv) TYPE(coeff_type), INTENT(IN) :: v_rspace + TYPE(global_environment_type), INTENT(IN) :: globenv TYPE(real_matrix_set_type), INTENT(INOUT) :: h - INTEGER, OPTIONAL, INTENT(OUT) :: nproduct + INTEGER, INTENT(OUT) :: nproduct ! *** Local variables *** @@ -86,7 +89,7 @@ CONTAINS INTEGER, DIMENSION(:), POINTER :: orb_neighbors REAL(wp), DIMENSION(:), POINTER :: orb_r2,rpgfa,rpgfb,zeta,zetb - REAL(wp), DIMENSION(:,:), POINTER :: h_block,hab,orb_r,work + REAL(wp), DIMENSION(:,:), POINTER :: h_block,hab,orb_r,s_block,work ! --------------------------------------------------------------------------- @@ -144,6 +147,8 @@ CONTAINS block_col=jatom,& block=h_block) + IF (.NOT.ASSOCIATED(h_block)) CYCLE + DO iset=1,nseta radius_set_a = basis_set_a%set_radius(iset) @@ -236,10 +241,14 @@ CONTAINS END DO - IF (PRESENT(nproduct)) nproduct = npgf_product +! *** Release work storage *** DEALLOCATE (hab,work) + CALL mp_sum(npgf_product,globenv%group) + + nproduct = npgf_product + CALL timestop(0.0_wp,handle) END SUBROUTINE integrate_v_rspace diff --git a/src/lib/fast.F b/src/lib/fast.F index 9d41c51c15..45a95e567a 100644 --- a/src/lib/fast.F +++ b/src/lib/fast.F @@ -24,13 +24,13 @@ SUBROUTINE rankup ( n, za, cmat, zb, ex, ey, ez, scr ) n3 = n2 * n ( 3 ) scr ( 1:n2 ) = CMPLX ( 0._dbl, KIND = dbl ) #if defined (__sp_lib) - CALL CGERU ( n ( 1 ), n ( 2 ), zb, ex, 1, ey, 1, scr, n ( 1 ) ) - CALL CSCAL ( n3, za, cmat, 1 ) - CALL CGERU ( n2, n ( 3 ), cone, scr, 1, ez, 1, cmat, n2 ) + CALL cgeru ( n ( 1 ), n ( 2 ), zb, ex, 1, ey, 1, scr, n ( 1 ) ) + CALL cscal ( n3, za, cmat, 1 ) + CALL cgeru ( n2, n ( 3 ), cone, scr, 1, ez, 1, cmat, n2 ) #else - CALL ZGERU ( n ( 1 ), n ( 2 ), zb, ex, 1, ey, 1, scr, n ( 1 ) ) - CALL ZSCAL ( n3, za, cmat, 1 ) - CALL ZGERU ( n2, n ( 3 ), cone, scr, 1, ez, 1, cmat, n2 ) + CALL zgeru ( n ( 1 ), n ( 2 ), zb, ex, 1, ey, 1, scr, n ( 1 ) ) + CALL zscal ( n3, za, cmat, 1 ) + CALL zgeru ( n2, n ( 3 ), cone, scr, 1, ez, 1, cmat, n2 ) #endif END SUBROUTINE rankup diff --git a/src/lib/mltfftsg.F b/src/lib/mltfftsg.F index cad4e920cc..68640590cf 100644 --- a/src/lib/mltfftsg.F +++ b/src/lib/mltfftsg.F @@ -34,47 +34,47 @@ SUBROUTINE mltfftsg ( transa, transb, a, ldax, lday, b, ldbx, ldby, & ISIG = -ISIGN TSCAL = ( ABS ( SCALE -1._dbl ) > 1.e-12_dbl ) - CALL CTRIG ( N, TRIG, AFTER, BEFORE, NOW, ISIG, IC ) + CALL ctrig ( N, TRIG, AFTER, BEFORE, NOW, ISIG, IC ) LOT = NCACHE / ( 4 * N ) LOT = LOT - MOD ( LOT + 1, 2 ) LOT = MAX ( 1, LOT ) DO ITR = 1, M, LOT NFFT = MIN ( M - ITR + 1, LOT ) IF ( TRANSA == 'N' .OR. TRANSA == 'n' ) THEN - CALL FFTPRE ( NFFT, NFFT, LDAX, LOT, N, A ( 1, ITR ), Z ( 1, 1 ), & + CALL fftpre ( NFFT, NFFT, LDAX, LOT, N, A ( 1, ITR ), Z ( 1, 1 ), & TRIG, NOW ( 1 ), AFTER ( 1 ), BEFORE ( 1 ), ISIG ) ELSE - CALL FFTSTP ( LDAX, NFFT, N, LOT, N, A ( ITR, 1 ), Z ( 1, 1 ), & + CALL fftstp ( LDAX, NFFT, N, LOT, N, A ( ITR, 1 ), Z ( 1, 1 ), & TRIG, NOW ( 1 ), AFTER ( 1 ), BEFORE ( 1 ), ISIG ) ENDIF IF ( TSCAL ) THEN IF ( LOT == NFFT ) THEN - CALL DSCAL ( 2 * LOT * N, SCALE, Z ( 1, 1 ), 1 ) + CALL dscal ( 2 * LOT * N, SCALE, Z ( 1, 1 ), 1 ) ELSE DO I = 1, N - CALL DSCAL ( 2 * NFFT, SCALE, Z ( LOT * ( I - 1 ) + 1, 1 ), 1 ) + CALL dscal ( 2 * NFFT, SCALE, Z ( LOT * ( I - 1 ) + 1, 1 ), 1 ) END DO END IF END IF IF(IC.EQ.1) THEN IF(TRANSB == 'N'.OR.TRANSB == 'n') THEN - CALL ZGETMO(Z(1,1),LOT,NFFT,N,B(1,ITR),LDBX) + CALL zgetmo(Z(1,1),LOT,NFFT,N,B(1,ITR),LDBX) ELSE - CALL MATMOV(NFFT,N,Z(1,1),LOT,B(ITR,1),LDBX) + CALL matmov(NFFT,N,Z(1,1),LOT,B(ITR,1),LDBX) ENDIF ELSE INZEE=1 DO I=2,IC-1 - CALL FFTSTP(LOT,NFFT,N,LOT,N,Z(1,INZEE), & + CALL fftstp(LOT,NFFT,N,LOT,N,Z(1,INZEE), & Z(1,3-INZEE),TRIG,NOW(I),AFTER(I), & BEFORE(I),ISIG) INZEE=3-INZEE ENDDO IF(TRANSB == 'N'.OR.TRANSB == 'n') THEN - CALL FFTROT(LOT,NFFT,N,NFFT,LDBX,Z(1,INZEE), & + CALL fftrot(LOT,NFFT,N,NFFT,LDBX,Z(1,INZEE), & B(1,ITR),TRIG,NOW(IC),AFTER(IC),BEFORE(IC),ISIG) ELSE - CALL FFTSTP(LOT,NFFT,N,LDBX,N,Z(1,INZEE), & + CALL fftstp(LOT,NFFT,N,LDBX,N,Z(1,INZEE), & B(ITR,1),TRIG,NOW(IC),AFTER(IC),BEFORE(IC),ISIG) ENDIF ENDIF diff --git a/src/library_tests.F b/src/library_tests.F index c56c48ef41..a49fc15734 100644 --- a/src/library_tests.F +++ b/src/library_tests.F @@ -381,7 +381,7 @@ SUBROUTINE matmul_test ( globenv ) tstart = m_cputime ( ) DO j = 1, ntim - CALL DGEMM ( "N", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len ) + CALL dgemm ( "N", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len ) END DO tend = m_cputime ( ) t = tend - tstart @@ -394,8 +394,8 @@ SUBROUTINE matmul_test ( globenv ) tstart = m_cputime ( ) DO j = 1, ntim - CALL DGEMM ( "N", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len ) - CALL DCOPY ( len * len , mc, 1, ma, 1) + CALL dgemm ( "N", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len ) + CALL dcopy ( len * len , mc, 1, ma, 1) END DO tend = m_cputime ( ) t = tend - tstart @@ -408,7 +408,7 @@ SUBROUTINE matmul_test ( globenv ) tstart = m_cputime ( ) DO j = 1, ntim - CALL DGEMM ( "N", "T", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len ) + CALL dgemm ( "N", "T", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len ) END DO tend = m_cputime ( ) t = tend - tstart @@ -421,7 +421,7 @@ SUBROUTINE matmul_test ( globenv ) tstart = m_cputime ( ) DO j = 1, ntim - CALL DGEMM ( "T", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len ) + CALL dgemm ( "T", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len ) END DO tend = m_cputime ( ) t = tend - tstart diff --git a/src/mathlib.F b/src/mathlib.F index e128e4f141..33de3406af 100644 --- a/src/mathlib.F +++ b/src/mathlib.F @@ -568,7 +568,7 @@ CONTAINS ! *** Externals (LAPACK) *** - REAL(wp), EXTERNAL :: dlange +!MK REAL(wp), EXTERNAL :: dlange EXTERNAL dgecon,dgerfs,dgetrf,dgetrs @@ -648,7 +648,7 @@ CONTAINS norm = 'I' END IF - a_norm = dlange(norm,n,n,a(:,:),n,work(:)) +!MK a_norm = dlange(norm,n,n,a(:,:),n,work(:)) ! *** Compute the reciprocal of the condition number of a *** diff --git a/src/matrix_types.F b/src/matrix_types.F index 380ecc8e0d..fd459644ce 100644 --- a/src/matrix_types.F +++ b/src/matrix_types.F @@ -31,17 +31,21 @@ MODULE matrix_types ! matrix_name,matrix_symmetry) ! SUBROUTINE allocate_full_real_matrix(matrix,nrow,ncol,matrix_name,& ! matrix_symmetry) +! SUBROUTINE copy_sparse_to_full_matrix(sparse_matrix,full_matrix) ! SUBROUTINE copy_real_matrix(source,target) ! SUBROUTINE deallocate_real_matrix(matrix) ! SUBROUTINE deallocate_real_matrix_row(matrix,block_row) -! SUBROUTINE get_matrix_info(matrix,nrow,ncol,nblock_row,nblock_col,& -! matrix_name,matrix_symmetry) +! SUBROUTINE get_matrix_info(matrix,matrix_name,matrix_symmetry,& +! nblock_row,nblock_col,nrow,ncol,& +! first_row,last_row,first_col,last_col,& +! nblock,nelement) ! SUBROUTINE get_real_block_node(block_node,block_col,block) ! SUBROUTINE get_real_matrix_block(matrix,block_row,block_col,block_node,block) ! SUBROUTINE put_real_block_node(block_node,matrix,block_row,block_col,block) ! SUBROUTINE put_real_matrix_block(matrix,block_row,block_col,block) ! SUBROUTINE replicate_real_matrix(source,target,target_name) ! SUBROUTINE replicate_real_matrix_structure(source,target,target_name) +! SUBROUTINE symmetrise_diagonal_blocks(matrix,option) ! FUNCTION find_real_block_node(matrix,block_row,block_col) RESULT(block_node) ! FUNCTION first_real_block_node(matrix,block_row) RESULT(first_block_node) @@ -51,6 +55,10 @@ MODULE matrix_types USE kinds, ONLY: wp => dp + USE memory_utilities, ONLY: reallocate + USE termination, ONLY: stop_memory,stop_program + USE timings, ONLY: timeset,timestop + IMPLICIT NONE PRIVATE @@ -98,13 +106,15 @@ MODULE matrix_types PUBLIC :: add_block_node,& allocate_matrix,& copy_matrix,& + copy_sparse_to_full_matrix,& deallocate_matrix,& deallocate_matrix_row,& get_block_node,& get_matrix_info,& put_block_node,& replicate_matrix,& - replicate_matrix_structure + replicate_matrix_structure,& + symmetrise_diagonal_blocks ! *** Public functions *** @@ -175,9 +185,6 @@ CONTAINS ! *************************************************************************** - USE memory_utilities, ONLY: reallocate - USE termination, ONLY: stop_memory,stop_program - TYPE(real_matrix_type), POINTER :: matrix INTEGER, INTENT(IN) :: block_col,block_row REAL(wp), DIMENSION(:,:), OPTIONAL, POINTER :: block @@ -243,9 +250,6 @@ CONTAINS ! *************************************************************************** - USE memory_utilities, ONLY: reallocate - USE termination, ONLY: stop_memory - TYPE(real_matrix_type), POINTER :: matrix CHARACTER(LEN=*), INTENT(IN) :: matrix_name,matrix_symmetry INTEGER, INTENT(IN) :: nblock_row,nblock_col,ncol,nrow @@ -321,9 +325,6 @@ CONTAINS ! *************************************************************************** - USE memory_utilities, ONLY: reallocate - USE termination, ONLY: stop_memory - TYPE(real_matrix_type), POINTER :: matrix CHARACTER(LEN=*), INTENT(IN) :: matrix_name,matrix_symmetry INTEGER, INTENT(IN) :: ncol,nrow @@ -384,6 +385,86 @@ CONTAINS END SUBROUTINE allocate_full_real_matrix +! ***************************************************************************** + + SUBROUTINE copy_sparse_to_full_matrix(sparse_matrix,full_matrix) + +! Purpose: Copy the matrix blocks of a sparse matrix to the corresponding +! full matrix which is allocated in this routine. + +! History: - Creation (19.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(real_matrix_type), POINTER :: sparse_matrix + REAL(wp), DIMENSION(:,:), POINTER :: full_matrix + +! *** Local parameters *** + + CHARACTER(LEN=*), PARAMETER :: routine =& + "SUBROUTINE copy_sparse_to_full_matrix (MODULE matrix_types)" + +! *** Local variables *** + + TYPE(real_block_node_type), POINTER :: block_node + + INTEGER :: iblock_col,iblock_row,icol,icol_block,irow,irow_block + + REAL(wp), DIMENSION(:,:), POINTER :: sparse_block + +! --------------------------------------------------------------------------- + +! *** Check the association status of the input matrix *** + + IF (.NOT.ASSOCIATED(sparse_matrix)) THEN + CALL stop_program(routine,"The input matrix pointer is not associated") + END IF + + IF (ASSOCIATED(full_matrix)) DEALLOCATE (full_matrix) + + full_matrix => reallocate(full_matrix,1,sparse_matrix%nrow,& + 1,sparse_matrix%ncol) + +! *** Traverse all block nodes of the sparse matrix *** + + DO iblock_row=1,sparse_matrix%nblock_row + + block_node => first_block_node(sparse_matrix,iblock_row) + + DO WHILE (ASSOCIATED(block_node)) + + CALL get_block_node(block_node=block_node,& + block_col=iblock_col,& + block=sparse_block) + + icol_block = 1 + + DO icol=sparse_matrix%first_col(iblock_col),& + sparse_matrix%last_col(iblock_col) + + irow_block = 1 + + DO irow=sparse_matrix%first_row(iblock_row),& + sparse_matrix%last_row(iblock_row) + + full_matrix(irow,icol) = sparse_block(irow_block,icol_block) + + irow_block = irow_block + 1 + + END DO + + icol_block = icol_block + 1 + + END DO + + block_node => next_block_node(block_node) + + END DO + + END DO + + END SUBROUTINE copy_sparse_to_full_matrix + ! ***************************************************************************** SUBROUTINE copy_real_matrix(source,target) @@ -394,8 +475,6 @@ CONTAINS ! *************************************************************************** - USE termination, ONLY: stop_program - TYPE(real_matrix_type), POINTER :: source,target ! *** Local parameters *** @@ -408,7 +487,8 @@ CONTAINS TYPE(real_block_node_type), POINTER :: source_block_node,& target_block_node - INTEGER :: first_col,first_row,iblock_col,iblock_row,icol,irow,istat,& + INTEGER :: first_col,first_row,handle,& + iblock_col,iblock_row,icol,irow,istat,& jblock_col,jblock_row,jcol,jrow,last_col,last_row,& source_first_col,source_first_row,& source_last_col,source_last_row,& @@ -423,6 +503,8 @@ CONTAINS CALL stop_program(routine,"The source matrix pointer is not associated") END IF + CALL timeset("copy_real_matrix","I","",handle) + ! *** Check the association status of the target matrix *** IF (ASSOCIATED(target)) THEN @@ -506,6 +588,8 @@ CONTAINS END IF + CALL timestop(0.0_wp,handle) + END SUBROUTINE copy_real_matrix ! ***************************************************************************** @@ -518,8 +602,6 @@ CONTAINS ! *************************************************************************** - USE termination, ONLY: stop_memory - TYPE(real_matrix_type), POINTER :: matrix ! *** Local parameters *** @@ -567,8 +649,6 @@ CONTAINS ! *************************************************************************** - USE termination, ONLY: stop_memory - TYPE(real_matrix_type), POINTER :: matrix INTEGER, INTENT(IN) :: block_row @@ -649,8 +729,10 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE get_matrix_info(matrix,nrow,ncol,nblock_row,nblock_col,& - matrix_name,matrix_symmetry) + SUBROUTINE get_matrix_info(matrix,matrix_name,matrix_symmetry,& + nblock_row,nblock_col,nrow,ncol,& + first_row,last_row,first_col,last_col,& + nblock_allocated,nelement_allocated) ! Purpose: Return the requested matrix information. @@ -658,12 +740,22 @@ CONTAINS ! *************************************************************************** - USE memory_utilities, ONLY: reallocate - TYPE(real_matrix_type), POINTER :: matrix CHARACTER(LEN=60), OPTIONAL, INTENT(OUT) :: matrix_name CHARACTER(LEN=40), OPTIONAL, INTENT(OUT) :: matrix_symmetry - INTEGER, OPTIONAL, INTENT(OUT) :: nblock_row,nblock_col,ncol,nrow + INTEGER, OPTIONAL, INTENT(OUT) :: nblock_allocated,nblock_row,& + nblock_col,ncol,& + nelement_allocated,nrow + INTEGER, DIMENSION(:), OPTIONAL, POINTER :: first_col,first_row,& + last_col,last_row + +! *** Local variables *** + + TYPE(real_block_node_type), POINTER :: block_node + + INTEGER :: iblock_row + + REAL(wp), DIMENSION(:,:), POINTER :: block ! --------------------------------------------------------------------------- @@ -673,6 +765,33 @@ CONTAINS IF (PRESENT(nblock_col)) nblock_col = matrix%nblock_col IF (PRESENT(nrow)) nrow = matrix%nrow IF (PRESENT(ncol)) ncol = matrix%ncol + IF (PRESENT(first_row)) first_row => matrix%first_row + IF (PRESENT(last_row)) last_row => matrix%last_row + IF (PRESENT(first_col)) first_col => matrix%first_col + IF (PRESENT(last_col)) last_col => matrix%last_col + + IF (PRESENT(nblock_allocated)) THEN + nblock_allocated = 0 + DO iblock_row=1,matrix%nblock_row + block_node => first_block_node(matrix,iblock_row) + DO WHILE (ASSOCIATED(block_node)) + nblock_allocated = nblock_allocated + 1 + block_node => next_block_node(block_node) + END DO + END DO + END IF + + IF (PRESENT(nelement_allocated)) THEN + nelement_allocated = 0 + DO iblock_row=1,matrix%nblock_row + block_node => first_block_node(matrix,iblock_row) + DO WHILE (ASSOCIATED(block_node)) + CALL get_block_node(block_node=block_node,block=block) + nelement_allocated = nelement_allocated + SIZE(block) + block_node => next_block_node(block_node) + END DO + END DO + END IF END SUBROUTINE get_matrix_info @@ -699,7 +818,9 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE get_real_matrix_block(matrix,block_row,block_col,block_node,block) + SUBROUTINE get_real_matrix_block(matrix,block_row,block_col,& + first_row,last_row,first_col,last_col,& + block_node,block) ! Purpose: Get block node data set. @@ -710,6 +831,8 @@ CONTAINS TYPE(real_matrix_type), POINTER :: matrix TYPE(real_block_node_type), OPTIONAL, POINTER :: block_node INTEGER, INTENT(IN) :: block_col,block_row + INTEGER, OPTIONAL, INTENT(OUT) :: first_col,first_row,& + last_col,last_row REAL(wp), DIMENSION(:,:), OPTIONAL, POINTER :: block ! *** Local variables *** @@ -718,6 +841,12 @@ CONTAINS ! --------------------------------------------------------------------------- + IF (PRESENT(first_row)) first_row = matrix%first_row(block_row) + IF (PRESENT(last_row)) last_row = matrix%last_row(block_row) + + IF (PRESENT(first_col)) first_col = matrix%first_col(block_col) + IF (PRESENT(last_col)) last_col = matrix%last_col(block_col) + current_block_node => find_real_block_node(matrix,block_row,block_col) IF (ASSOCIATED(current_block_node)) THEN @@ -760,9 +889,6 @@ CONTAINS ! *************************************************************************** - USE memory_utilities, ONLY: reallocate - USE termination, ONLY: stop_program - TYPE(real_block_node_type), POINTER :: block_node TYPE(real_matrix_type), POINTER :: matrix INTEGER, INTENT(IN) :: block_row @@ -814,8 +940,6 @@ CONTAINS ! *************************************************************************** - USE termination, ONLY: stop_program - TYPE(real_matrix_type), POINTER :: matrix INTEGER, INTENT(IN) :: block_col,block_row REAL(wp), DIMENSION(:,:), OPTIONAL, POINTER:: block @@ -887,8 +1011,6 @@ CONTAINS ! *************************************************************************** - USE termination, ONLY: stop_program - TYPE(real_matrix_type), POINTER :: source,target CHARACTER(LEN=*), INTENT(IN) :: target_name @@ -913,6 +1035,8 @@ CONTAINS CALL stop_program(routine,"The source matrix pointer is not associated") END IF + IF (ASSOCIATED(target)) CALL deallocate_real_matrix(target) + ! *** Allocate a new matrix structure *** CALL allocate_real_matrix(matrix=target,& @@ -964,8 +1088,6 @@ CONTAINS ! *************************************************************************** - USE termination, ONLY: stop_program - TYPE(real_matrix_type), POINTER :: source,target CHARACTER(LEN=*), INTENT(IN) :: target_name @@ -1025,6 +1147,67 @@ CONTAINS END SUBROUTINE replicate_real_matrix_structure +! ***************************************************************************** + + SUBROUTINE symmetrise_diagonal_blocks(matrix) + +! Purpose: Symmetrise the diagonal blocks of matrix. + +! History: - Creation (13.06.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(real_matrix_type), POINTER :: matrix + +! *** Local parameters *** + + CHARACTER(LEN=*), PARAMETER :: routine =& + "SUBROUTINE symmetrise_diagonal_blocks (MODULE matrix_types)" + +! *** Local variables *** + + TYPE(real_block_node_type), POINTER :: block_node + + INTEGER :: iblock_col,iblock_row,icol,irow + + REAL(wp), DIMENSION(:,:), POINTER :: block + +! --------------------------------------------------------------------------- + +! *** Check the association status of the input matrix *** + + IF (.NOT.ASSOCIATED(matrix)) THEN + CALL stop_program(routine,"The input matrix pointer is not associated") + END IF + +! *** Traverse all block nodes *** + + DO iblock_row=1,matrix%nblock_row + + block_node => first_block_node(matrix,iblock_row) + + DO WHILE (ASSOCIATED(block_node)) + + CALL get_block_node(block_node=block_node,& + block_col=iblock_col,& + block=block) + + IF (iblock_row == iblock_col) THEN + DO irow=1,SIZE(block,1) + DO icol=irow+1,SIZE(block,2) + block(icol,irow) = block(irow,icol) + END DO + END DO + END IF + + block_node => next_block_node(block_node) + + END DO + + END DO + + END SUBROUTINE symmetrise_diagonal_blocks + ! ***************************************************************************** END MODULE matrix_types diff --git a/src/message_passing.F b/src/message_passing.F index 823b058034..54c81b5f02 100644 --- a/src/message_passing.F +++ b/src/message_passing.F @@ -46,7 +46,7 @@ MODULE message_passing INTERFACE mp_sum MODULE PROCEDURE mp_sum_i1, mp_sum_r1, mp_sum_c1, mp_sum_iv, & mp_sum_rv, mp_sum_cv, mp_sum_im, mp_sum_rm, mp_sum_cm, & - mp_sum_im3, mp_sum_rm3, mp_sum_root_rm3 + mp_sum_im3, mp_sum_rm3, mp_sum_cm3, mp_sum_root_rm, mp_sum_root_rm3 END INTERFACE INTERFACE mp_max @@ -1079,6 +1079,34 @@ SUBROUTINE mp_sum_rm(msg,gid) mp_perf ( 3 ) % time = mp_perf ( 3 ) % time + ( t_end - t_start ) #endif END SUBROUTINE mp_sum_rm +SUBROUTINE mp_sum_root_rm(msg,root,gid) + IMPLICIT NONE + REAL ( dbl ), INTENT ( INOUT ) :: msg ( :, : ) + INTEGER, INTENT ( IN ) :: root,gid + INTEGER :: msglen, m1, m2, ierr, taskid + REAL ( dbl ), ALLOCATABLE :: res ( :, : ) +#if defined(__parallel) + t_start = cputime ( ) + CALL mpi_comm_rank ( gid, taskid, ierr ) + IF ( ierr /= 0 ) CALL mp_stop( ierr, "mpi_comm_rank @ mp_sum_root_rm" ) + msglen = SIZE(msg) + m1 = SIZE(msg,1) + m2 = SIZE(msg,2) + ALLOCATE (res(m1,m2),STAT=ierr) + IF ( ierr /= 0 ) CALL mp_stop( 0, "allocate @ mp_sum_root_rm" ) + CALL mpi_reduce(msg,res,msglen,MPI_DOUBLE_PRECISION,MPI_SUM,& + root,gid,ierr) + IF ( taskid == root ) THEN + msg = res + END IF + DEALLOCATE (res) + IF ( ierr /= 0 ) CALL mp_stop( ierr, "mpi_reduce @ mp_sum_root_rm" ) + mp_perf ( 3 ) % count = mp_perf ( 3 ) % count + 1 + mp_perf ( 3 ) % msg_size = mp_perf ( 3 ) % msg_size + msglen * reallen + t_end = cputime ( ) + mp_perf ( 3 ) % time = mp_perf ( 3 ) % time + ( t_end - t_start ) +#endif +END SUBROUTINE mp_sum_root_rm SUBROUTINE mp_sum_rm3(msg,gid) IMPLICIT NONE REAL ( dbl ), INTENT ( INOUT ) :: msg ( :, :, : ) @@ -1126,7 +1154,7 @@ SUBROUTINE mp_sum_root_rm3(msg,root,gid) msg = res END IF DEALLOCATE (res) - IF ( ierr /= 0 ) CALL mp_stop( ierr, "mpi_reduce @ mp_sum_root_rm" ) + IF ( ierr /= 0 ) CALL mp_stop( ierr, "mpi_reduce @ mp_sum_root_rm3" ) mp_perf ( 3 ) % count = mp_perf ( 3 ) % count + 1 mp_perf ( 3 ) % msg_size = mp_perf ( 3 ) % msg_size + msglen * reallen t_end = cputime ( ) @@ -1199,6 +1227,31 @@ SUBROUTINE mp_sum_cm(msg,gid) mp_perf ( 3 ) % time = mp_perf ( 3 ) % time + ( t_end - t_start ) #endif END SUBROUTINE mp_sum_cm +SUBROUTINE mp_sum_cm3(msg,gid) + IMPLICIT NONE + COMPLEX ( dbl ), INTENT ( INOUT ) :: msg ( :, :, : ) + INTEGER, INTENT ( IN ) :: gid + INTEGER :: msglen, m1, m2, m3, ierr + COMPLEX ( dbl ), ALLOCATABLE :: res ( :, :, : ) +#if defined(__parallel) + t_start = cputime ( ) + msglen = 2*size(msg) + m1 = size(msg,1) + m2 = size(msg,2) + m3 = size(msg,3) + ALLOCATE (res(m1,m2,m3),STAT=ierr) + IF ( ierr /= 0 ) CALL mp_stop( 0, "allocate @ mp_sum_cm3" ) + CALL mpi_allreduce(msg,res,msglen,MPI_DOUBLE_PRECISION,MPI_SUM,gid, & + ierr) + msg = res + DEALLOCATE (res) + IF ( ierr /= 0 ) CALL mp_stop( ierr, "mpi_allreduce @ mp_sum_cm3" ) + mp_perf ( 3 ) % count = mp_perf ( 3 ) % count + 1 + mp_perf ( 3 ) % msg_size = mp_perf ( 3 ) % msg_size + msglen * reallen*2 + t_end = cputime ( ) + mp_perf ( 3 ) % time = mp_perf ( 3 ) % time + ( t_end - t_start ) +#endif +END SUBROUTINE mp_sum_cm3 !****************************************************************************** diff --git a/src/mo_types.F b/src/mo_types.F index 962b5956b5..bbe6ceefa6 100644 --- a/src/mo_types.F +++ b/src/mo_types.F @@ -25,13 +25,20 @@ MODULE mo_types USE kinds, ONLY: wp => dp + USE blacs, ONLY: allocate_blacs_matrix,& + blacs_matrix_type,& + copy_blacs_to_full_matrix,& + deallocate_blacs_matrix,& + get_blacs_info,& + read_blacs_matrix,& + write_blacs_matrix USE global_types, ONLY: global_environment_type USE input_utilities, ONLY: close_file,& open_file USE matrix_types, ONLY: allocate_matrix,& - deallocate_matrix,& get_block_node,& - real_matrix_set_type + real_matrix_type + USE message_passing, ONLY: mp_sum IMPLICIT NONE @@ -40,7 +47,7 @@ MODULE mo_types TYPE mo_set_type INTEGER :: homo,lfomo REAL(wp), DIMENSION(:), POINTER :: eigenvalues,occupation_numbers - TYPE(real_matrix_set_type) :: eigenvectors + TYPE(blacs_matrix_type) :: eigenvectors END TYPE mo_set_type ! *** Public data types *** @@ -66,12 +73,21 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE allocate_mo_set(mo_set,nao,nmo) + SUBROUTINE allocate_mo_set(mo_set,nao,nmo,nrow_block,ncol_block,context,& + globenv) + +! Purpose: Allocate a wavefunction data structure. + +! History: - Creation (15.05.2001, Matthias Krack) + +! *************************************************************************** USE memory_utilities, ONLY: reallocate - TYPE(mo_set_type), INTENT(OUT) :: mo_set - INTEGER, INTENT(IN) :: nao,nmo + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(mo_set_type), INTENT(OUT) :: mo_set + INTEGER, INTENT(IN) :: context,nao,ncol_block,nmo,& + nrow_block ! --------------------------------------------------------------------------- @@ -84,12 +100,14 @@ CONTAINS NULLIFY (mo_set%occupation_numbers) mo_set%occupation_numbers => reallocate(mo_set%occupation_numbers,1,nmo) - NULLIFY (mo_set%eigenvectors%matrix) - CALL allocate_matrix(matrix=mo_set%eigenvectors%matrix,& - nrow=nao,& - ncol=nmo,& - matrix_name="MO EIGENVECTORS",& - matrix_symmetry="none") + CALL allocate_blacs_matrix(new_matrix=mo_set%eigenvectors,& + nrow_global=nao,& + ncol_global=nmo,& + nrow_block=nrow_block,& + ncol_block=ncol_block,& + name="MO EIGENVECTORS",& + context=context,& + globenv=globenv) END SUBROUTINE allocate_mo_set @@ -102,42 +120,56 @@ CONTAINS ! --------------------------------------------------------------------------- DEALLOCATE (mo_set%eigenvalues) - DEALLOCATE (mo_set%occupation_numbers) - - CALL deallocate_matrix(matrix=mo_set%eigenvectors%matrix) + CALL deallocate_blacs_matrix(mo_set%eigenvectors) END SUBROUTINE deallocate_mo_set ! ***************************************************************************** - SUBROUTINE read_mo_set(mo_set,globenv) + SUBROUTINE read_mo_set(mo_set,context,globenv) + +! Purpose: Read the MO eigenvectors from the restart file. + +! History: - Creation (15.05.2001, Matthias Krack) +! - Parallel input (19.05.2001, MK) + +! *************************************************************************** TYPE(global_environment_type), INTENT(IN) :: globenv - TYPE(mo_set_type), INTENT(IN) :: mo_set + TYPE(mo_set_type), INTENT(OUT) :: mo_set + INTEGER, INTENT(IN) :: context ! *** Local variables *** - INTEGER :: iao,imo,restart_unit - - REAL(wp), DIMENSION(:,:), POINTER :: mo_eigenvectors + CHARACTER(LEN=6) :: extension + CHARACTER(LEN=200) :: file_name + INTEGER :: mype,npe,restart_unit ! --------------------------------------------------------------------------- - CALL get_block_node(matrix=mo_set%eigenvectors%matrix,& - block_row=1,& - block_col=1,& - block=mo_eigenvectors) + CALL get_blacs_info(context=context,& + globenv=globenv,& + my_process_number=mype,& + number_of_processes=npe) - CALL open_file(file_name=globenv%restart_file_name,& + IF (npe > 1) THEN + WRITE (UNIT=extension,FMT="(I6)") mype + file_name = TRIM(globenv%restart_file_name)//"."//ADJUSTL(extension) + ELSE + file_name = globenv%restart_file_name + END IF + + CALL open_file(file_name=file_name,& file_action="READ",& file_form="UNFORMATTED",& file_status="OLD",& unit_number=restart_unit) - READ (UNIT=restart_unit) ((mo_eigenvectors(iao,imo),& - iao=1,SIZE(mo_eigenvectors,1)),& - imo=1,SIZE(mo_eigenvectors,2)) + CALL read_blacs_matrix(matrix=mo_set%eigenvectors,& + lunit=restart_unit,& + context=context,& + globenv=globenv) CALL close_file(unit_number=restart_unit) @@ -145,33 +177,48 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE write_mo_set_to_restart_unit(mo_set,globenv) + SUBROUTINE write_mo_set_to_restart_unit(mo_set,context,globenv) + +! Purpose: Write the MO eigenvectors to the restart file. + +! History: - Creation (15.05.2001, Matthias Krack) + +! *************************************************************************** TYPE(global_environment_type), INTENT(IN) :: globenv TYPE(mo_set_type), INTENT(IN) :: mo_set + INTEGER, INTENT(IN) :: context ! *** Local variables *** - INTEGER :: iao,imo,restart_unit - - REAL(wp), DIMENSION(:,:), POINTER :: mo_eigenvectors + CHARACTER(LEN=6) :: extension + CHARACTER(LEN=200) :: file_name + INTEGER :: mype,npe,restart_unit ! --------------------------------------------------------------------------- - CALL get_block_node(matrix=mo_set%eigenvectors%matrix,& - block_row=1,& - block_col=1,& - block=mo_eigenvectors) + CALL get_blacs_info(context=context,& + globenv=globenv,& + my_process_number=mype,& + number_of_processes=npe) - CALL open_file(file_name=globenv%restart_file_name,& + IF (npe > 1) THEN + WRITE (UNIT=extension,FMT="(I6)") mype + file_name = TRIM(globenv%restart_file_name)//"."//ADJUSTL(extension) + ELSE + file_name = globenv%restart_file_name + END IF + + CALL open_file(file_name=file_name,& file_action="WRITE",& file_form="UNFORMATTED",& file_status="REPLACE",& unit_number=restart_unit) - WRITE (UNIT=restart_unit) ((mo_eigenvectors(iao,imo),& - iao=1,SIZE(mo_eigenvectors,1)),& - imo=1,SIZE(mo_eigenvectors,2)) + CALL write_blacs_matrix(matrix=mo_set%eigenvectors,& + lunit=restart_unit,& + context=context,& + globenv=globenv) CALL close_file(unit_number=restart_unit) @@ -179,7 +226,7 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE write_mo_set_to_output_unit(mo_set,before,after,globenv) + SUBROUTINE write_mo_set_to_output_unit(mo_set,before,after,context,globenv) ! Purpose: Write the MO eigenvalues, MO occupation numbers and ! MO eigenvectors. @@ -200,8 +247,8 @@ CONTAINS USE orbital_pointers, ONLY: nso TYPE(global_environment_type), INTENT(IN) :: globenv - TYPE(mo_set_type), INTENT(IN) :: mo_set - INTEGER, INTENT(IN) :: after,before + TYPE(mo_set_type), INTENT(INOUT) :: mo_set + INTEGER, INTENT(IN) :: after,before,context ! *** Local variables *** @@ -220,10 +267,11 @@ CONTAINS output_unit = globenv%scr - CALL get_block_node(matrix=mo_set%eigenvectors%matrix,& - block_row=1,& - block_col=1,& - block=matrix) + IF (globenv%print%mo_eigenvectors) THEN + NULLIFY (matrix) + CALL copy_blacs_to_full_matrix(mo_set%eigenvectors,matrix,context,& + globenv) + END IF IF (.NOT.globenv%ionode) RETURN @@ -317,6 +365,10 @@ CONTAINS WRITE (output_unit,"(/)") +! *** Release work storage *** + + DEALLOCATE (matrix) + END SUBROUTINE write_mo_set_to_output_unit ! ***************************************************************************** diff --git a/src/om_utilities.F b/src/om_utilities.F index 6c86faa9d6..00ffe0f5f1 100644 --- a/src/om_utilities.F +++ b/src/om_utilities.F @@ -37,7 +37,8 @@ MODULE om_utilities ! *** Public subroutines *** - PUBLIC :: write_cartesian_matrix,& + PUBLIC :: write_blacs_matrix,& + write_cartesian_matrix,& write_spherical_matrix,& write_g_matrix @@ -47,7 +48,7 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE write_cartesian_matrix(block_matrix,before,after,globenv) + SUBROUTINE write_cartesian_matrix(sparse_matrix,before,after,globenv) ! Purpose: Write a Cartesian matrix. @@ -70,29 +71,23 @@ CONTAINS USE atoms, ONLY: atom_info,natom,ncgf USE global_types, ONLY: global_environment_type USE mathlib, ONLY: symmetrize_matrix - USE matrix_types, ONLY: add_block_node,& - allocate_matrix,& - copy_matrix,& - deallocate_matrix,& - get_block_node,& + USE matrix_types, ONLY: copy_sparse_to_full_matrix,& get_matrix_info,& real_matrix_type USE orbital_pointers, ONLY: nco TYPE(global_environment_type), INTENT(IN) :: globenv - TYPE(real_matrix_type), POINTER :: block_matrix + TYPE(real_matrix_type), POINTER :: sparse_matrix INTEGER, INTENT(IN) :: after,before ! *** Local variables *** - TYPE(real_matrix_type), POINTER :: full_matrix - CHARACTER(LEN=60) :: name - CHARACTER(LEN=40) :: symmetry - CHARACTER(LEN=25) :: fmtstr1 - CHARACTER(LEN=35) :: fmtstr2 - INTEGER :: from,iatom,icgf,ico,icol,ikind,irow,& - iset,ishell,jcol,l,left,ncol,& - output_unit,right,to,width + CHARACTER(LEN=60) :: name + CHARACTER(LEN=40) :: symmetry + CHARACTER(LEN=25) :: fmtstr1 + CHARACTER(LEN=35) :: fmtstr2 + INTEGER :: from,iatom,icgf,ico,icol,ikind,irow,iset,ishell,jcol,& + l,left,ncol,output_unit,right,to,width REAL(wp), DIMENSION(:,:), POINTER :: matrix @@ -100,24 +95,12 @@ CONTAINS output_unit = globenv%scr - CALL get_matrix_info(matrix=block_matrix,& + CALL get_matrix_info(matrix=sparse_matrix,& matrix_name=name,& matrix_symmetry=symmetry) - NULLIFY (full_matrix) - - CALL allocate_matrix(matrix=full_matrix,& - nrow=ncgf,& - ncol=ncgf,& - matrix_name="WORK",& - matrix_symmetry=symmetry) - - CALL copy_matrix(block_matrix,full_matrix) - - CALL get_block_node(matrix=full_matrix,& - block_row=1,& - block_col=1,& - block=matrix) + NULLIFY (matrix) + CALL copy_sparse_to_full_matrix(sparse_matrix,matrix) IF (symmetry == "symmetric") THEN CALL symmetrize_matrix(matrix,"upper_to_lower") @@ -125,9 +108,7 @@ CONTAINS CALL symmetrize_matrix(matrix,"anti_upper_to_lower") END IF -#if defined(__parallel) CALL mp_sum(matrix,globenv%group) -#endif IF (.NOT.globenv%ionode) RETURN @@ -187,13 +168,13 @@ CONTAINS ! *** Release work storage *** - CALL deallocate_matrix(full_matrix) + DEALLOCATE (matrix) END SUBROUTINE write_cartesian_matrix ! ***************************************************************************** - SUBROUTINE write_spherical_matrix(block_matrix,before,after,globenv,& + SUBROUTINE write_spherical_matrix(sparse_matrix,before,after,globenv,& matrix_name) ! Purpose: Write a spherical matrix. @@ -217,30 +198,24 @@ CONTAINS USE atoms, ONLY: atom_info,natom,nsgf USE global_types, ONLY: global_environment_type USE mathlib, ONLY: symmetrize_matrix - USE matrix_types, ONLY: add_block_node,& - allocate_matrix,& - deallocate_matrix,& - copy_matrix,& - get_block_node,& + USE matrix_types, ONLY: copy_sparse_to_full_matrix,& get_matrix_info,& real_matrix_type USE orbital_pointers, ONLY: nso TYPE(global_environment_type), INTENT(IN) :: globenv - TYPE(real_matrix_type), POINTER :: block_matrix + TYPE(real_matrix_type), POINTER :: sparse_matrix INTEGER, INTENT(IN) :: after,before CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: matrix_name ! *** Local variables *** - TYPE(real_matrix_type), POINTER :: full_matrix - CHARACTER(LEN=60) :: name - CHARACTER(LEN=40) :: symmetry - CHARACTER(LEN=25) :: fmtstr1 - CHARACTER(LEN=35) :: fmtstr2 - INTEGER :: from,iatom,icol,ikind,irow,iset,isgf,& - iso,ishell,jcol,l,left,ncol,& - output_unit,right,to,width + CHARACTER(LEN=60) :: name + CHARACTER(LEN=40) :: symmetry + CHARACTER(LEN=25) :: fmtstr1 + CHARACTER(LEN=35) :: fmtstr2 + INTEGER :: from,iatom,icol,ikind,irow,iset,isgf,iso,ishell,jcol,& + l,left,ncol,output_unit,right,to,width REAL(wp), DIMENSION(:,:), POINTER :: matrix @@ -248,26 +223,14 @@ CONTAINS output_unit = globenv%scr - CALL get_matrix_info(matrix=block_matrix,& + CALL get_matrix_info(matrix=sparse_matrix,& matrix_name=name,& matrix_symmetry=symmetry) IF (PRESENT(matrix_name)) name = matrix_name - NULLIFY (full_matrix) - - CALL allocate_matrix(matrix=full_matrix,& - nrow=nsgf,& - ncol=nsgf,& - matrix_name="WORK",& - matrix_symmetry=symmetry) - - CALL copy_matrix(block_matrix,full_matrix) - - CALL get_block_node(matrix=full_matrix,& - block_row=1,& - block_col=1,& - block=matrix) + NULLIFY (matrix) + CALL copy_sparse_to_full_matrix(sparse_matrix,matrix) IF (symmetry == "symmetric") THEN CALL symmetrize_matrix(matrix,"upper_to_lower") @@ -275,9 +238,7 @@ CONTAINS CALL symmetrize_matrix(matrix,"anti_upper_to_lower") END IF -#if defined(__parallel) CALL mp_sum(matrix,globenv%group) -#endif IF (.NOT.globenv%ionode) RETURN @@ -337,13 +298,13 @@ CONTAINS ! *** Release work storage *** - CALL deallocate_matrix(full_matrix) + DEALLOCATE (matrix) END SUBROUTINE write_spherical_matrix ! ***************************************************************************** - SUBROUTINE write_g_matrix(block_matrix,before,after,globenv) + SUBROUTINE write_g_matrix(sparse_matrix,before,after,globenv) ! Purpose: Write the Cartesian G matrix. @@ -366,29 +327,23 @@ CONTAINS USE atoms, ONLY: atom_info,natom,ncgf_aux USE global_types, ONLY: global_environment_type USE mathlib, ONLY: symmetrize_matrix - USE matrix_types, ONLY: add_block_node,& - allocate_matrix,& - copy_matrix,& - deallocate_matrix,& - get_block_node,& + USE matrix_types, ONLY: copy_sparse_to_full_matrix,& get_matrix_info,& real_matrix_type USE orbital_pointers, ONLY: nco TYPE(global_environment_type), INTENT(IN) :: globenv - TYPE(real_matrix_type), POINTER :: block_matrix + TYPE(real_matrix_type), POINTER :: sparse_matrix INTEGER, INTENT(IN) :: after,before ! *** Local variables *** - TYPE(real_matrix_type), POINTER :: full_matrix - CHARACTER(LEN=60) :: name - CHARACTER(LEN=40) :: symmetry - CHARACTER(LEN=25) :: fmtstr1 - CHARACTER(LEN=35) :: fmtstr2 - INTEGER :: from,iatom,icgf,ico,icol,ikind,irow,& - iset,ishell,jcol,l,left,ncol,& - output_unit,right,to,width + CHARACTER(LEN=60) :: name + CHARACTER(LEN=40) :: symmetry + CHARACTER(LEN=25) :: fmtstr1 + CHARACTER(LEN=35) :: fmtstr2 + INTEGER :: from,iatom,icgf,ico,icol,ikind,irow,iset,ishell,jcol,& + l,left,ncol,output_unit,right,to,width REAL(wp), DIMENSION(:,:), POINTER :: matrix @@ -396,24 +351,13 @@ CONTAINS output_unit = globenv%scr - CALL get_matrix_info(matrix=block_matrix,& + CALL get_matrix_info(matrix=sparse_matrix,& matrix_name=name,& matrix_symmetry=symmetry) - NULLIFY (full_matrix) - CALL allocate_matrix(matrix=full_matrix,& - nrow=ncgf_aux,& - ncol=ncgf_aux,& - matrix_name="WORK",& - matrix_symmetry=symmetry) - - CALL copy_matrix(block_matrix,full_matrix) - - CALL get_block_node(matrix=full_matrix,& - block_row=1,& - block_col=1,& - block=matrix) + NULLIFY (matrix) + CALL copy_sparse_to_full_matrix(sparse_matrix,matrix) IF (symmetry == "symmetric") THEN CALL symmetrize_matrix(matrix,"upper_to_lower") @@ -421,9 +365,7 @@ CONTAINS CALL symmetrize_matrix(matrix,"anti_upper_to_lower") END IF -#if defined(__parallel) CALL mp_sum(matrix,globenv%group) -#endif IF (.NOT.globenv%ionode) RETURN @@ -483,10 +425,120 @@ CONTAINS ! *** Release work storage *** - CALL deallocate_matrix(full_matrix) + DEALLOCATE (matrix) END SUBROUTINE write_g_matrix +! ***************************************************************************** + + SUBROUTINE write_blacs_matrix(blacs_matrix,before,after,context,globenv,& + matrix_name) + +! Purpose: Write a spherical matrix. + +! History: - Creation (12.06.2001, Matthias Krack) + +! *************************************************************************** + + USE atomic_kinds, ONLY: kind_info + USE atoms, ONLY: atom_info,natom,nsgf + USE blacs, ONLY: blacs_matrix_type,& + copy_blacs_to_full_matrix,& + get_blacs_matrix_info + USE global_types, ONLY: global_environment_type + USE orbital_pointers, ONLY: nso + + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(blacs_matrix_type), INTENT(INOUT) :: blacs_matrix + INTEGER, INTENT(IN) :: after,before,context + CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: matrix_name + +! *** Local variables *** + + CHARACTER(LEN=60) :: name + CHARACTER(LEN=25) :: fmtstr1 + CHARACTER(LEN=35) :: fmtstr2 + INTEGER :: from,iatom,icol,ikind,irow,iset,isgf,iso,ishell,jcol,& + l,left,ncol,ncol_global,nrow_global,output_unit,& + right,to,width + + REAL(wp), DIMENSION(:,:), POINTER :: matrix + +! --------------------------------------------------------------------------- + + output_unit = globenv%scr + + NULLIFY (matrix) + CALL copy_blacs_to_full_matrix(blacs_matrix,matrix,context,globenv) + + IF (.NOT.globenv%ionode) RETURN + + IF (PRESENT(matrix_name)) THEN + name = matrix_name + ELSE + CALL get_blacs_matrix_info(blacs_matrix,name=name) + END IF + +! *** Definition of the variable formats *** + + fmtstr1 = "(/,T2,23X, ( X,I5, X))" + fmtstr2 = "(T2,2I5,2X,A2,1X,A8, (1X,F . ))" + +! *** Write headline *** + + WRITE (output_unit,"(/,/,T2,A)") TRIM(name) + +! *** Write the variable format strings *** + + width = before + after + 3 + ncol = INT(56/width) + + right = MAX((after-2),1) + left = width - right - 5 + + WRITE (fmtstr1(11:12),"(I2)") ncol + WRITE (fmtstr1(14:15),"(I2)") left + WRITE (fmtstr1(21:22),"(I2)") right + + WRITE (fmtstr2(22:23),"(I2)") ncol + WRITE (fmtstr2(29:30),"(I2)") width - 1 + WRITE (fmtstr2(32:33),"(I2)") after + +! *** Write the matrix in the selected format *** + + DO icol=1,nsgf,ncol + from = icol + to = MIN((from+ncol-1),nsgf) + WRITE (output_unit,fmtstr1) (jcol,jcol=from,to) + irow = 1 + DO iatom=1,natom + IF (iatom /= 1) WRITE (output_unit,"(A)") "" + ikind = atom_info(iatom)%kind + isgf = 1 + DO iset=1,kind_info(ikind)%orb_basis_set%nset + DO ishell=1,kind_info(ikind)%orb_basis_set%nshell(iset) + l = kind_info(ikind)%orb_basis_set%l(ishell,iset) + DO iso=1,nso(l) + WRITE (output_unit,fmtstr2)& + irow,iatom,kind_info(ikind)%element_symbol,& + kind_info(ikind)%orb_basis_set%sgf_symbol(isgf),& + (matrix(irow,jcol),jcol=from,to) + isgf = isgf + 1 + irow = irow + 1 + END DO + END DO + END DO + END DO + END DO + + WRITE (output_unit,"(/)") + +! *** Release work storage *** + + DEALLOCATE (matrix) + + END SUBROUTINE write_blacs_matrix + ! ***************************************************************************** END MODULE om_utilities diff --git a/src/print_keys.F b/src/print_keys.F index aaed287470..877ae0cfd9 100644 --- a/src/print_keys.F +++ b/src/print_keys.F @@ -32,6 +32,7 @@ MODULE print_keys atomic_coordinates,& basic_data_types,& basis_set,& + blacs_info,& cartesian_om,& cell_parameters,& charge_density_matrix,& @@ -133,6 +134,7 @@ CONTAINS print_key%atomic_coordinates = key_value(level,(/t,t,t/)) print_key%basic_data_types = key_value(level,(/f,f,f/)) print_key%basis_set = key_value(level,(/f,t,t/)) + print_key%blacs_info = key_value(level,(/f,t,t/)) print_key%cartesian_om = key_value(level,(/f,f,f/)) print_key%cell_parameters = key_value(level,(/t,t,t/)) print_key%charge_density_matrix = key_value(level,(/f,f,t/)) diff --git a/src/qs_scf.F b/src/qs_scf.F index ba4820ffd4..030cc43b89 100644 --- a/src/qs_scf.F +++ b/src/qs_scf.F @@ -28,9 +28,28 @@ MODULE qs_scf USE atomic_kinds, ONLY: kind_info,nkind USE atoms, ONLY: atom_info USE basis_set_types, ONLY: gto_basis_set_type,maxlcgf + USE blacs, ONLY: allocate_blacs_matrix,& + blacs_gemm,& + blacs_matrix_type,& + blacs_set_all,& + blacs_set_element,& + blacs_syevx,& + blacs_symm,& + blacs_syrk,& + copy_blacs_to_sparse_matrix,& + copy_sparse_to_blacs_matrix,& + finish_blacs,& + get_blacs_matrix_info,& + power_blacs_matrix,& + replicate_blacs_matrix,& + start_blacs,& + symmetrise_blacs_matrix USE coefficient_types, ONLY: coeff_allocate,& coeff_type,& coeff_zero + USE core_energies, ONLY: calculate_ecore,& + calculate_ecore_overlap,& + calculate_ecore_self USE core_hamiltonian, ONLY: build_core_hamiltonian_matrix,h,s USE diis, ONLY: eps_diis,max_diis,scf_diis USE functionals, ONLY: init_functionals,vwn_c,vwn_x @@ -39,24 +58,25 @@ MODULE qs_scf read_object,& search,& start_input_session - USE mathlib, ONLY: diagonalize_matrix,& - power_matrix,& - symmetrize_matrix USE matrix_types, ONLY: allocate_matrix,& copy_matrix,& + deallocate_matrix,& first_block_node,& get_block_node,& get_matrix_info,& next_block_node,& real_block_node_type,& real_matrix_set_type,& - replicate_matrix_structure + replicate_matrix_structure,& + symmetrise_diagonal_blocks USE memory_utilities, ONLY: reallocate + USE message_passing, ONLY: mp_max,mp_sum,mp_sync USE mo_types, ONLY: allocate_mo_set,& mo_set_type,& read_mo_set,& write_mo_set - USE om_utilities, ONLY: write_spherical_matrix + USE om_utilities, ONLY: write_blacs_matrix,& + write_spherical_matrix USE pw_grid_types, ONLY: HALFSPACE,pw_grid_type USE pw_grids, ONLY: pw_find_cutoff,& pw_grid_construct,& @@ -73,8 +93,9 @@ MODULE qs_scf PRIVATE TYPE(coeff_type) :: rho_gspace,rho_rspace,v_rspace + TYPE(blacs_matrix_type) :: ortho,scf_work1,scf_work2 TYPE(mo_set_type) :: alpha_mo - TYPE(real_matrix_set_type) :: ks,ortho,p,scf_work1,scf_work2 + TYPE(real_matrix_set_type) :: ks,p TYPE(pw_grid_type) :: pw_grid CHARACTER(LEN=10) :: density_guess = "CORE" @@ -83,6 +104,7 @@ MODULE qs_scf ecore_overlap = 0.0_wp,& ecore_self = 0.0_wp,& ehartree = 0.0_wp,& + eps_eigval = 1.0E-5_wp,& eps_scf = 1.0E-5_wp,& etotal = 0.0_wp,& ex = 0.0_wp,& @@ -90,9 +112,14 @@ MODULE qs_scf total_rho_core_rspace = 0.0_wp,& total_rho_elec_rspace = 0.0_wp,& total_rho_gspace = 0.0_wp,& - total_rho_rspace = 0.0_wp + total_rho_rspace = 0.0_wp,& + work_syevx = 0.0_wp INTEGER :: max_scf = 50,& - nelectron = 0 + nelectron = 0,& + nrow_block = 32,& + ncol_block = 32,& + nprow = 0,& + npcol = 0 ! *** Public variables *** @@ -122,42 +149,36 @@ CONTAINS ! *** Local variables *** REAL(wp) :: delta,diis_error,t1,t2 - INTEGER :: handle,iscf,output_unit - LOGICAL :: diis_step - - REAL(wp), DIMENSION(:,:), POINTER :: scf_work1_matrix + INTEGER :: context,handle,iscf,output_unit + LOGICAL :: diis_step,ionode ! --------------------------------------------------------------------------- - CALL timeset("scf","I","",handle) - - output_unit = globenv%scr - ! *** Quick return, if no SCF iteration is requested *** IF (max_scf < 1) RETURN + CALL timeset("scf","I","",handle) + + ionode = globenv%ionode + output_unit = globenv%scr + CALL init_functionals() CALL init_grid(globenv) - IF (globenv%print%scf) THEN + IF (ionode.AND.globenv%print%scf) THEN WRITE (UNIT=output_unit,FMT="(/,/,T2,A)")& "SCF WAVEFUNCTION OPTIMIZATION" END IF - CALL init_scf_run(globenv,alpha_mo) - - CALL get_block_node(matrix=scf_work1%matrix,& - block_row=1,& - block_col=1,& - block=scf_work1_matrix) + CALL init_scf_run(alpha_mo,context,globenv) etotal = 0.0_wp iscf = 0 diis_step = .FALSE. - IF (globenv%print%scf) THEN + IF (ionode.AND.globenv%print%scf) THEN WRITE (UNIT=output_unit,& FMT="(/,T3,A,T9,A,T34,A,T49,A,T68,A,/,T3,A)")& "Step","Update method","Time","Convergence","Total energy",& @@ -180,39 +201,40 @@ CONTAINS CALL write_spherical_matrix(ks%matrix,4,6,globenv) END IF - scf_work1_matrix(:,:) = 0.0_wp - CALL copy_matrix(ks%matrix,scf_work1%matrix) - CALL symmetrize_matrix(scf_work1_matrix,"upper_to_lower") + CALL copy_sparse_to_blacs_matrix(ks%matrix,scf_work1,context,globenv) +!MK if only gemm is used +! CALL symmetrise_blacs_matrix(scf_work1,scf_work2,context,globenv) IF (iscf > 1) THEN - CALL scf_diis(globenv,alpha_mo,scf_work1,scf_work2,delta,diis_error,& - diis_step) + CALL scf_diis(alpha_mo,scf_work1,scf_work2,delta,diis_error,diis_step,& + context,globenv) END IF - CALL orthogonalize_matrix(ortho,scf_work1,scf_work2) + CALL orthogonalise_matrix(ortho,scf_work1,scf_work2,context,globenv) - CALL eigensolver(scf_work1,alpha_mo,ortho) + CALL eigensolver(scf_work1,alpha_mo,ortho,scf_work2,context,globenv) IF (globenv%print%each_scf_step) THEN - CALL write_mo_set(alpha_mo,4,6,globenv) + CALL write_mo_set(alpha_mo,4,6,context,globenv) END IF - CALL calculate_density_matrix(alpha_mo,scf_work1) + CALL calculate_density_matrix(alpha_mo,scf_work1,context,globenv) - CALL copy_matrix(scf_work1%matrix,ks%matrix) + CALL copy_blacs_to_sparse_matrix(scf_work1,ks%matrix,context,globenv) + CALL symmetrise_diagonal_blocks(ks%matrix) t2 = cputime() IF (diis_step) THEN - CALL density_mixing(ks,p,1.0_wp,delta) - IF (globenv%print%scf) THEN + CALL density_mixing(ks,p,1.0_wp,delta,globenv) + IF (ionode.AND.globenv%print%scf) THEN WRITE (UNIT=output_unit,& FMT="(T2,I5,2X,A,T15,E10.2,T30,F8.2,T40,2F20.10)")& iscf,"DIIS",diis_error,t2 - t1,delta,etotal END IF ELSE - CALL density_mixing(ks,p,p_mix,delta) - IF (globenv%print%scf) THEN + CALL density_mixing(ks,p,p_mix,delta,globenv) + IF (ionode.AND.globenv%print%scf) THEN WRITE (UNIT=output_unit,& FMT="(T2,I5,2X,A,T15,F6.2,T30,F8.2,T40,2F20.10)")& iscf,"Mixing",p_mix,t2 - t1,delta,etotal @@ -220,14 +242,13 @@ CONTAINS END IF IF (delta < eps_scf) THEN - IF (globenv%print%scf) THEN + IF (ionode.AND.globenv%print%scf) THEN WRITE(UNIT=output_unit,FMT="(/,T3,A,/)")& "*** SCF run converged ***" END IF - IF (.NOT.diis_step) CALL copy_matrix(ks%matrix,p%matrix) EXIT ELSE IF (iscf == max_scf) THEN - IF (globenv%print%scf) THEN + IF (ionode.AND.globenv%print%scf) THEN WRITE(UNIT=output_unit,FMT="(/,T3,A,/)")& "*** SCF run NOT converged ***" END IF @@ -236,7 +257,7 @@ CONTAINS END DO - IF (globenv%print%scf) THEN + IF (ionode.AND.globenv%print%scf) THEN WRITE (UNIT=output_unit,FMT="(/,(T3,A,T40,2F20.10))")& "Total electronic density (r-space): ",& total_rho_elec_rspace,total_rho_elec_rspace + REAL(nelectron,wp),& @@ -258,7 +279,7 @@ CONTAINS IF (.NOT.diis_step) CALL copy_matrix(ks%matrix,p%matrix) - CALL write_mo_set(alpha_mo,4,6,globenv) + CALL write_mo_set(alpha_mo,4,6,context,globenv) IF (globenv%print%density_matrix) THEN CALL write_spherical_matrix(p%matrix,4,6,globenv) @@ -271,7 +292,9 @@ CONTAINS ! *** Write restart file *** - CALL write_mo_set(alpha_mo,globenv) + CALL write_mo_set(alpha_mo,context,globenv) + + CALL finish_blacs(context,globenv) CALL timestop(0.0_wp,handle) @@ -284,7 +307,6 @@ CONTAINS USE collocate_density, ONLY: calculate_rho_core,& calculate_rho_elec,& calculate_total_rho - USE core_energies, ONLY: calculate_ecore USE integrate_potential, ONLY: integrate_v_rspace USE print_keys, ONLY: DEBUG @@ -293,29 +315,32 @@ CONTAINS ! *** Local variables *** INTEGER :: handle,nproduct,output_unit + LOGICAL :: ionode ! --------------------------------------------------------------------------- CALL timeset("build_kohn_sham_matrix","I","",handle) + ionode = globenv%ionode output_unit = globenv%scr - CALL calculate_ecore(h,p,ecore) + CALL calculate_ecore(h,p,ecore,globenv) CALL coeff_zero(rho_rspace) - CALL calculate_rho_elec(p,rho_rspace,total_rho_elec_rspace,nproduct) + CALL calculate_rho_elec(p,rho_rspace,total_rho_elec_rspace,nproduct,& + globenv) - IF (globenv%print%total_densities) THEN + IF (ionode.AND.globenv%print%total_densities) THEN WRITE (UNIT=output_unit,FMT="(/,T3,A,I10)")& "Number of collocated products (r-space):",nproduct END IF CALL coeff_zero(v_rspace) - CALL calculate_xc_potential(rho_rspace%pw,v_rspace%pw) + CALL calculate_xc_potential(rho_rspace,v_rspace,globenv) - CALL calculate_rho_core(rho_rspace,total_rho_rspace) + CALL calculate_rho_core(rho_rspace,total_rho_rspace,globenv) total_rho_core_rspace = total_rho_rspace - total_rho_elec_rspace @@ -324,7 +349,7 @@ CONTAINS total_rho_gspace = calculate_total_rho(rho_gspace) - IF (globenv%print%total_densities) THEN + IF (ionode.AND.globenv%print%total_densities) THEN WRITE (UNIT=output_unit,FMT="(/,(T3,A,T40,2F20.10))")& "Total electronic density (r-space): ",& total_rho_elec_rspace,total_rho_elec_rspace + REAL(nelectron,wp),& @@ -335,11 +360,11 @@ CONTAINS "Total charge density (g-space): ",total_rho_gspace END IF - CALL calculate_hartree_potential(rho_gspace) + CALL calculate_hartree_potential(rho_gspace,globenv) etotal = ecore_overlap + ecore_self + ecore + ehartree + ex + ec - IF (globenv%print%scf_energies) THEN + IF (ionode.AND.globenv%print%scf_energies) THEN WRITE (UNIT=output_unit,FMT="(/,(T3,A,T60,F20.10))")& "Overlap energy of the core charge distribution:",ecore_overlap,& "Self energy of the core charge distribution: ",ecore_self,& @@ -357,9 +382,9 @@ CONTAINS CALL copy_matrix(h%matrix,ks%matrix) - CALL integrate_v_rspace(v_rspace,ks,nproduct) + CALL integrate_v_rspace(v_rspace,ks,nproduct,globenv) - IF (globenv%print%total_densities) THEN + IF (ionode.AND.globenv%print%total_densities) THEN WRITE (UNIT=output_unit,FMT="(/,T3,A,I10)")& "Number of integrated products (r-space):",nproduct END IF @@ -379,22 +404,26 @@ CONTAINS ! *** Local variables *** INTEGER :: handle,output_unit + LOGICAL :: ionode ! --------------------------------------------------------------------------- CALL timeset("init_grid","I","",handle) + ionode = globenv%ionode output_unit = globenv%scr CALL pw_grid_construct(pw_grid) pw_grid%grid_span = HALFSPACE - IF (globenv%print%pw_grid_information) THEN - CALL pw_grid_setup(cell,pw_grid,cutoff,info=output_unit,& + IF (ionode.AND.globenv%print%pw_grid_information) THEN + CALL pw_grid_setup(cell,pw_grid,cutoff,& + info=output_unit,& orthorhombic=.TRUE.) ELSE - CALL pw_grid_setup(cell,pw_grid,cutoff,orthorhombic=.TRUE.) + CALL pw_grid_setup(cell,pw_grid,cutoff,& + orthorhombic=.TRUE.) END IF CALL coeff_allocate(rho_rspace,pw_grid,REALDATA3D) @@ -412,19 +441,17 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE init_scf_run(globenv,mo) + SUBROUTINE init_scf_run(mo,context,globenv) -! Purpose: Initialize a SCF run. +! Purpose: Initialise a SCF run. ! History: - Creation (30.04.2001, Matthias Krack) ! *************************************************************************** - USE core_energies, ONLY: calculate_ecore_overlap,& - calculate_ecore_self - TYPE(global_environment_type), INTENT(IN) :: globenv TYPE(mo_set_type), INTENT(OUT) :: mo + INTEGER, INTENT(OUT) :: context ! *** Local parameters *** @@ -433,15 +460,24 @@ CONTAINS ! *** Local variables *** - CHARACTER(LEN=40) :: symmetry - INTEGER :: handle,homo,ikind,nao,ndep,nmo,output_unit + INTEGER :: handle,homo,ikind,nao,ndep,nmo,output_unit + LOGICAL :: ionode + + INTEGER, DIMENSION(0:globenv%num_pe-1) :: nblock_pe,nelement_pe ! --------------------------------------------------------------------------- CALL timeset("init_scf_run","I","",handle) + ionode = globenv%ionode output_unit = globenv%scr +! *** Initialise BLACS *** + + CALL start_blacs(nprow,npcol,context,globenv) + +! *** Calculate the number of electrons *** + nelectron = 0 DO ikind=1,nkind @@ -457,58 +493,76 @@ CONTAINS ! *** Get the dimension of the full SCF matrices, *** ! *** i.e. the total number of atomic orbitals *** - CALL get_matrix_info(matrix=h%matrix,& - nrow=nao,& - matrix_symmetry=symmetry) + CALL get_matrix_info(matrix=h%matrix,nrow=nao) - nmo = nao -!MK nmo = homo + IF (globenv%print%mo_eigenvectors) THEN + nmo = nao + ELSE + nmo = homo + END IF - CALL allocate_mo_set(mo,nao,nmo) +! *** Allocate the distributed MO eigenvectors *** + + CALL allocate_mo_set(mo,nao,nmo,nrow_block,ncol_block,context,globenv) mo%homo = homo mo%occupation_numbers(1:homo) = 2.0_wp -! *** Allocate the full SCF matrices *** +! *** Get BLACS block size of the MO eigenvector matrix *** +! *** which has to fit to the other distributed SCF matrices *** - NULLIFY (ortho%matrix) - CALL allocate_matrix(matrix=ortho%matrix,& - nrow=nao,& - ncol=nao,& - matrix_name="ORTHOGONALIZATION MATRIX",& - matrix_symmetry=symmetry) + CALL get_blacs_matrix_info(matrix=mo%eigenvectors,& + nrow_block=nrow_block,& + ncol_block=ncol_block) - NULLIFY (scf_work1%matrix) - CALL allocate_matrix(matrix=scf_work1%matrix,& - nrow=nao,& - ncol=nao,& - matrix_name="SCF WORK MATRIX 1",& - matrix_symmetry="none") +! *** Allocate the distributed SCF matrices *** - NULLIFY (scf_work2%matrix) - CALL allocate_matrix(matrix=scf_work2%matrix,& - nrow=nao,& - ncol=nao,& - matrix_name="SCF WORK MATRIX 2",& - matrix_symmetry="none") + CALL allocate_blacs_matrix(new_matrix=ortho,& + nrow_global=nao,& + ncol_global=nao,& + nrow_block=nrow_block,& + ncol_block=ncol_block,& + name="ORTHOGONALIZATION MATRIX",& + context=context,& + globenv=globenv) + + CALL replicate_blacs_matrix(prototype_matrix=ortho,& + new_matrix=scf_work1,& + name="SCF WORK MATRIX 1") + + CALL replicate_blacs_matrix(prototype_matrix=ortho,& + new_matrix=scf_work2,& + name="SCF WORK MATRIX 2") + + CALL copy_sparse_to_blacs_matrix(h%matrix,scf_work1,context,globenv) + +! *** Redistribute the core Hamiltonian matrix *** + +!MK IF (gpw) THEN + CALL deallocate_matrix(h%matrix) + CALL replicate_matrix_structure(s%matrix,h%matrix,& + "CORE HAMILTONIAN MATRIX") + CALL copy_blacs_to_sparse_matrix(scf_work1,h%matrix,context,globenv) +!MK END IF NULLIFY (ks%matrix) CALL replicate_matrix_structure(h%matrix,ks%matrix,"KOHN-SHAM MATRIX") + NULLIFY (p%matrix) CALL replicate_matrix_structure(h%matrix,p%matrix,"DENSITY MATRIX") CALL calculate_ecore_self(ecore_self) - CALL calculate_ecore_overlap(ecore_overlap) + CALL calculate_ecore_overlap(ecore_overlap,globenv) - IF (globenv%print%scf_energies) THEN + IF (ionode.AND.globenv%print%scf_energies) THEN WRITE (UNIT=output_unit,FMT="(/,(T3,A,T60,F20.10))")& "Self energy of the core charge distribution: ",ecore_self,& "Overlap energy of the core charge distribution:",ecore_overlap END IF - CALL calculate_ortho_matrix(ortho,scf_work1,ndep) + CALL calculate_ortho_matrix(ortho,scf_work2,ndep,context,globenv) - IF (globenv%print%scf) THEN + IF (ionode.AND.globenv%print%scf) THEN WRITE (UNIT=output_unit,FMT="(/,(T3,A,I10))")& "Number of electrons: ",nelectron,& "Number of occupied orbitals: ",homo,& @@ -517,10 +571,11 @@ CONTAINS END IF IF (globenv%print%ortho_matrix) THEN - CALL write_spherical_matrix(ortho%matrix,4,6,globenv) + CALL write_blacs_matrix(ortho,4,6,context,globenv) END IF - CALL calculate_first_density_matrix(ortho,mo,p,scf_work1,scf_work2,globenv) + CALL calculate_first_density_matrix(ortho,mo,p,scf_work1,scf_work2,& + context,globenv) CALL timestop(0.0_wp,handle) @@ -528,43 +583,32 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE calculate_ortho_matrix(ortho,work,ndep) + SUBROUTINE calculate_ortho_matrix(ortho,work,ndep,context,globenv) - TYPE(real_matrix_set_type), INTENT(OUT) :: ortho - TYPE(real_matrix_set_type), INTENT(INOUT) :: work +! Purpose: Calculate the orthogonalization matrix (S**(-1/2)) + +! History: - Creation (01.05.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(blacs_matrix_type), INTENT(OUT) :: ortho + TYPE(blacs_matrix_type), INTENT(INOUT) :: work + INTEGER, INTENT(IN) :: context INTEGER, INTENT(OUT) :: ndep ! *** Local variables *** INTEGER :: handle,nao - REAL(wp), DIMENSION(:,:), POINTER :: ortho_matrix,work_matrix - ! --------------------------------------------------------------------------- CALL timeset("calculate_ortho_matrix","I","",handle) - CALL get_matrix_info(matrix=ortho%matrix,nrow=nao) - - CALL get_block_node(matrix=ortho%matrix,& - block_row=1,& - block_col=1,& - block=ortho_matrix) - - CALL get_block_node(matrix=work%matrix,& - block_row=1,& - block_col=1,& - block=work_matrix) - -! *** Calculate the orthogonalization matrix (S**(-1/2)) *** - - work_matrix(:,:) = 0.0_wp - - CALL copy_matrix(s%matrix,work%matrix) - - CALL symmetrize_matrix(work_matrix,"upper_to_lower") - - CALL power_matrix(work_matrix,ortho_matrix,-0.5_wp,1.0E-5_wp,ndep,.FALSE.) + CALL copy_sparse_to_blacs_matrix(s%matrix,ortho,context,globenv) + CALL power_blacs_matrix(ortho,work,-0.5_wp,eps_eigval,ndep,work_syevx,& + context,globenv) + CALL symmetrise_blacs_matrix(ortho,work,context,globenv) CALL timestop(0.0_wp,handle) @@ -572,92 +616,74 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE orthogonalize_matrix(ortho,ks,work) + SUBROUTINE orthogonalise_matrix(ortho,ks,work,context,globenv) - TYPE(real_matrix_set_type), INTENT(IN) :: ortho - TYPE(real_matrix_set_type), INTENT(INOUT) :: ks,work +! Purpose: Orthogonaliseation matrix (S**(-1/2)) + +! History: - Creation (01.05.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(IN) :: ortho + TYPE(blacs_matrix_type), INTENT(INOUT) :: ks,work + TYPE(global_environment_type), INTENT(IN) :: globenv + INTEGER, INTENT(IN) :: context ! *** Local variables *** INTEGER :: handle,nao - REAL(wp), DIMENSION(:,:), POINTER :: ks_matrix,ortho_matrix,work_matrix - ! --------------------------------------------------------------------------- - CALL timeset("orthogonalize_matrix","I","",handle) + CALL timeset("orthogonalise_matrix","I","",handle) - CALL get_matrix_info(matrix=ks%matrix,nrow=nao) - - CALL get_block_node(matrix=ks%matrix,& - block_row=1,& - block_col=1,& - block=ks_matrix) - - CALL get_block_node(matrix=ortho%matrix,& - block_row=1,& - block_col=1,& - block=ortho_matrix) - - CALL get_block_node(matrix=work%matrix,& - block_row=1,& - block_col=1,& - block=work_matrix) - - CALL dgemm("T","N",nao,nao,nao,1.0_wp,ks_matrix,nao,ortho_matrix,nao,& - 0.0_wp,work_matrix,nao) - - CALL dgemm("N","N",nao,nao,nao,1.0_wp,ortho_matrix,nao,work_matrix,nao,& - 0.0_wp,ks_matrix,nao) + CALL get_blacs_matrix_info(matrix=ks,nrow_global=nao) + CALL blacs_symm("L","U",nao,nao,1.0_wp,ks,ortho,0.0_wp,work,context,& + globenv) +!MK CALL blacs_gemm("T","N",nao,nao,nao,1.0_wp,ks,ortho,0.0_wp,work,context,& +!MK globenv) + CALL blacs_gemm("N","N",nao,nao,nao,1.0_wp,ortho,work,0.0_wp,ks,context,& + globenv) CALL timestop(0.0_wp,handle) - END SUBROUTINE orthogonalize_matrix + END SUBROUTINE orthogonalise_matrix ! ***************************************************************************** - SUBROUTINE eigensolver(ks,mo,ortho) + SUBROUTINE eigensolver(ks,mo,ortho,work,context,globenv) - TYPE(real_matrix_set_type), INTENT(IN) :: ortho - TYPE(real_matrix_set_type), INTENT(INOUT) :: ks +! Purpose: Diagonalise the Kohn-Sham matrix to get a new set of MO eigen- +! vectors and MO eigenvalues. + +! History: - Creation (01.05.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(IN) :: ortho + TYPE(blacs_matrix_type), INTENT(INOUT) :: ks,work + TYPE(global_environment_type), INTENT(IN) :: globenv TYPE(mo_set_type), INTENT(INOUT) :: mo + INTEGER, INTENT(IN) :: context ! *** Local variables *** INTEGER :: handle,nao,nmo - REAL(wp), DIMENSION(:,:), POINTER :: ks_matrix,mo_eigenvectors,ortho_matrix - ! --------------------------------------------------------------------------- CALL timeset("eigensolver","I","",handle) - CALL get_matrix_info(matrix=mo%eigenvectors%matrix,nrow=nao,ncol=nmo) + CALL get_blacs_matrix_info(matrix=mo%eigenvectors,& + nrow_global=nao,& + ncol_global=nmo) - CALL get_block_node(matrix=ks%matrix,& - block_row=1,& - block_col=1,& - block=ks_matrix) + CALL blacs_syevx(ks,work,mo%eigenvalues,nmo,work_syevx,context,globenv) - CALL get_block_node(matrix=mo%eigenvectors%matrix,& - block_row=1,& - block_col=1,& - block=mo_eigenvectors) - - CALL get_block_node(matrix=ortho%matrix,& - block_row=1,& - block_col=1,& - block=ortho_matrix) - -! CALL diagonalize_matrix(ks_matrix,mo_eigenvectors,mo%eigenvalues) -! CALL dgemm("N","N",nao,nmo,nao,1.0_wp,ortho_matrix,nao,mo_eigenvectors,& -! nao,0.0_wp,ks_matrix,nao) -! mo_eigenvectors(:,:) = ks_matrix(:,:) - - CALL diagonalize_matrix(ks_matrix,mo%eigenvalues,.FALSE.) - - CALL dgemm("N","N",nao,nmo,nao,1.0_wp,ortho_matrix,nao,ks_matrix,& - nao,0.0_wp,mo_eigenvectors,nao) +!MK CALL blacs_gemm("T","N",nao,nmo,nao,1.0_wp,ortho,work,0.0_wp,& +!MK mo%eigenvectors,context,globenv) + CALL blacs_symm("L","U",nao,nmo,1.0_wp,ortho,work,0.0_wp,mo%eigenvectors,& + context,globenv) CALL timestop(0.0_wp,handle) @@ -665,37 +691,30 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE calculate_density_matrix(mo,density) + SUBROUTINE calculate_density_matrix(mo,density,context,globenv) - TYPE(real_matrix_set_type), INTENT(OUT) :: density - TYPE(mo_set_type), INTENT(IN) :: mo +! Purpose: Calculate the density matrix from the MO eigenvectors and the +! MO occupation numbers. + +! History: - Creation (01.05.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(blacs_matrix_type), INTENT(OUT) :: density + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(mo_set_type), INTENT(IN) :: mo + INTEGER, INTENT(IN) :: context ! *** Local variables *** - INTEGER :: handle,nao - - REAL(wp), DIMENSION(:,:), POINTER :: mo_eigenvectors,density_matrix + INTEGER :: handle ! --------------------------------------------------------------------------- CALL timeset("calculate_density_matrix","I","",handle) - CALL get_matrix_info(matrix=mo%eigenvectors%matrix,nrow=nao) - - CALL get_block_node(matrix=mo%eigenvectors%matrix,& - block_row=1,& - block_col=1,& - block=mo_eigenvectors) - - CALL get_block_node(matrix=density%matrix,& - block_row=1,& - block_col=1,& - block=density_matrix) - - CALL dsyrk("U","N",nao,mo%homo,2.0_wp,mo_eigenvectors,nao,0.0_wp,& - density_matrix,nao) - - CALL symmetrize_matrix(density_matrix,"upper_to_lower") + CALL blacs_syrk("U","N",mo%homo,2.0_wp,mo%eigenvectors,0.0_wp,density,& + context,globenv) CALL timestop(0.0_wp,handle) @@ -703,8 +722,16 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE density_mixing(new_density,old_density,p_mix,delta) + SUBROUTINE density_mixing(new_density,old_density,p_mix,delta,globenv) +! Purpose: Perform a density mixing of the old (last SCF iteration) and the +! new (current) density matrix. + +! History: - Creation (01.05.2001, Matthias Krack) + +! *************************************************************************** + + TYPE(global_environment_type), INTENT(IN) :: globenv TYPE(real_matrix_set_type), INTENT(INOUT) :: new_density,old_density REAL(wp), INTENT(IN) :: p_mix REAL(wp), INTENT(OUT) :: delta @@ -757,17 +784,20 @@ CONTAINS END DO + CALL mp_max(delta,globenv%group) + CALL timestop(0.0_wp,handle) END SUBROUTINE density_mixing ! ***************************************************************************** - SUBROUTINE calculate_hartree_potential(rho) + SUBROUTINE calculate_hartree_potential(rho,globenv) USE mathconstants, ONLY: fourpi,twopi - TYPE(coeff_type), TARGET, INTENT(INOUT) :: rho + TYPE(coeff_type), TARGET, INTENT(INOUT) :: rho + TYPE(global_environment_type), INTENT(IN) :: globenv ! *** Local variables *** @@ -798,29 +828,28 @@ CONTAINS ehartree = 0.0_wp DO kg=lb_grid(3),ub_grid(3) - gz = REAL(kg,wp)*dg(3) IF (kg < 0) THEN k = kg + 1 + ng(3) ELSE k = kg + 1 END IF + gz = REAL(kg,wp)*dg(3) DO jg=lb_grid(2),ub_grid(2) - gy = REAL(jg,wp)*dg(2) IF (jg < 0) THEN j = jg + 1 + ng(2) ELSE j = jg + 1 END IF + gy = REAL(jg,wp)*dg(2) DO ig=lb_grid(1),ub_grid(1) - gx = REAL(ig,wp)*dg(1) IF (ig < 0) THEN i = ig + 1 + ng(1) ELSE i = ig + 1 END IF + gx = REAL(ig,wp)*dg(1) g2 = gx*gx + gy*gy + gz*gz IF (g2 > 1.0E-12_wp) THEN -!!!write (*,"(3I6,2F15.6)") i,j,k,grid(i,j,k) vhartree = fourpi*grid(i,j,k)/g2 ehartree = ehartree + REAL(grid(i,j,k))*REAL(vhartree) +& AIMAG(grid(i,j,k))*AIMAG(vhartree) @@ -838,27 +867,50 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE calculate_xc_potential(rho,vxc) + SUBROUTINE calculate_xc_potential(rho,vxc,globenv) - TYPE(pw_type), INTENT(IN) :: rho - TYPE(pw_type), INTENT(INOUT) :: vxc + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(coeff_type), TARGET, INTENT(IN) :: rho + TYPE(coeff_type), TARGET, INTENT(INOUT) :: vxc ! *** Local variables *** - INTEGER :: handle + INTEGER :: handle,k,kg + + INTEGER, DIMENSION(:), POINTER :: lb_grid,ng,ub_grid ! --------------------------------------------------------------------------- CALL timeset("calculate_xc_potential","I","",handle) + ng => rho%pw%pw_grid%npts(:) + lb_grid => rho%pw%pw_grid%bounds(1,:) + ub_grid => rho%pw%pw_grid%bounds(2,:) + ex = 0.0_wp ec = 0.0_wp - CALL vwn_x(rho%cr3d,ex,vxc%cr3d) - CALL vwn_c(rho%cr3d,ec,vxc%cr3d) + DO kg=lb_grid(3),ub_grid(3) - ex = rho%pw_grid%dvol*ex - ec = rho%pw_grid%dvol*ec + IF (kg < 0) THEN + k = kg + 1 + ng(3) + ELSE + k = kg + 1 + END IF + + IF (globenv%mepos /= MODULO(k,globenv%num_pe)) CYCLE + + CALL vwn_x(rho%pw%cr3d(:,:,kg),ex,vxc%pw%cr3d(:,:,kg)) + CALL vwn_c(rho%pw%cr3d(:,:,kg),ec,vxc%pw%cr3d(:,:,kg)) + + END DO + + CALL mp_sum(vxc%pw%cr3d,globenv%group) + CALL mp_sum(ex,globenv%group) + CALL mp_sum(ec,globenv%group) + + ex = rho%pw%pw_grid%dvol*ex + ec = rho%pw%pw_grid%dvol*ec CALL timestop(0.0_wp,handle) @@ -866,12 +918,15 @@ CONTAINS ! ***************************************************************************** - SUBROUTINE calculate_first_density_matrix(ortho,mo,p,work1,work2,globenv) + SUBROUTINE calculate_first_density_matrix(ortho,mo,p,work1,work2,context,& + globenv) + TYPE(blacs_matrix_type), INTENT(IN) :: ortho + TYPE(blacs_matrix_type), INTENT(INOUT) :: work1,work2 TYPE(global_environment_type), INTENT(IN) :: globenv - TYPE(real_matrix_set_type), INTENT(IN) :: ortho - TYPE(real_matrix_set_type), INTENT(INOUT) :: p,work1,work2 TYPE(mo_set_type), INTENT(INOUT) :: mo + TYPE(real_matrix_set_type), INTENT(INOUT) :: p + INTEGER, INTENT(IN) :: context ! *** Local variables *** @@ -883,38 +938,35 @@ CONTAINS INTEGER, DIMENSION(:), POINTER :: elec_conf - REAL(wp), DIMENSION(:,:), POINTER :: work1_matrix,work2_matrix - ! --------------------------------------------------------------------------- CALL timeset("calculate_first_density_matrix","I","",handle) -! *** Calculate an initial density guess (core Hamiltonian guess) *** - - CALL get_block_node(matrix=work1%matrix,& - block_row=1,& - block_col=1,& - block=work1_matrix) - - work1_matrix(:,:) = 0.0_wp + CALL blacs_set_all(work1,0.0_wp,context,globenv) IF (density_guess == "RESTART") THEN - CALL read_mo_set(mo,globenv) - CALL calculate_density_matrix(mo,work1) + CALL read_mo_set(mo,context,globenv) + CALL calculate_density_matrix(mo,work1,context,globenv) + CALL copy_blacs_to_sparse_matrix(work1,p%matrix,context,globenv) + CALL symmetrise_diagonal_blocks(p%matrix) ELSE IF (density_guess == "CORE") THEN - CALL copy_matrix(h%matrix,work1%matrix) - CALL symmetrize_matrix(work1_matrix,"upper_to_lower") - CALL orthogonalize_matrix(ortho,work1,work2) - CALL eigensolver(work1,mo,ortho) - CALL calculate_density_matrix(mo,work1) +! *** It is assumed that work1 holds the upper *** +! *** triangle part of the core Hamiltonian matrix *** + +!MK if only gemm is used +! CALL symmetrise_blacs_matrix(work1,work2,context,globenv) + CALL orthogonalise_matrix(ortho,work1,work2,context,globenv) + CALL eigensolver(work1,mo,ortho,work2,context,globenv) + CALL calculate_density_matrix(mo,work1,context,globenv) + CALL copy_blacs_to_sparse_matrix(work1,p%matrix,context,globenv) + CALL symmetrise_diagonal_blocks(p%matrix) ELSE IF (density_guess == "ATOMIC") THEN NULLIFY (elec_conf) - elec_conf => reallocate(elec_conf,0,maxlcgf) DO ikind=1,nkind @@ -949,7 +1001,7 @@ CONTAINS DO iatom=1,kind_info(ikind)%natom atoma = kind_info(ikind)%atom_list(iatom) isgf = atom_info(atoma)%first_sgf + isgfa - 1 - work1_matrix(isgf,isgf) = paa + CALL blacs_set_element(work1,isgf,isgf,paa,context,globenv) END DO END DO END IF @@ -958,14 +1010,17 @@ CONTAINS END DO + CALL copy_blacs_to_sparse_matrix(work1,p%matrix,context,globenv) + ELSE - STOP "qs_scf: wrong keyword for guess" + CALL stop_program("SUBROUTINE calculate_first_density_matrix "//& + "(MODULE qs_scf)",& + "An invalid keyword for the initial density "//& + "guess was specified") END IF - CALL copy_matrix(work1%matrix,p%matrix) - END SUBROUTINE calculate_first_density_matrix ! ***************************************************************************** @@ -996,11 +1051,18 @@ CONTAINS ! *** Load the default values *** density_guess = "CORE" + eps_eigval = 1.0E-5_wp eps_scf = 1.0E-5_wp eps_diis = 0.1_wp max_diis = 4 max_scf = 50 + nprow = 0 + npcol = 0 + nrow_block = 32 + ncol_block = 32 p_mix = 0.4_wp + work_syevx = 0.0_wp + globenv%restart_file_name = "RESTART" CALL start_input_session(globenv%input_file_name,globenv) @@ -1025,14 +1087,25 @@ CONTAINS CALL read_object(p_mix) CASE ("EPS_DIIS") CALL read_object(eps_diis) + CASE ("EPS_EIGVAL") + CALL read_object(eps_eigval) CASE ("EPS_SCF") CALL read_object(eps_scf) CASE ("MAX_DIIS") CALL read_object(max_diis) CASE ("MAX_SCF") CALL read_object(max_scf) + CASE ("BLOCKSIZE") + CALL read_object(nrow_block) + ncol_block = nrow_block + CASE ("PROCESS_GRID") + CALL read_object(nprow) + CALL read_object(npcol) CASE ("RESTART_FILE_NAME","RESTART_FILE","RESTART") CALL read_object(globenv%restart_file_name) + CASE ("WORK_SYEVX") + CALL read_object(work_syevx) + work_syevx = MIN(MAX(0.0_wp,work_syevx),1.0_wp) CASE DEFAULT IF (keyword == end_section) THEN EXIT @@ -1067,16 +1140,17 @@ CONTAINS IF (max_scf < 1) RETURN - WRITE (lunit,"(/,/,T2,A,/)") "SCF PARAMETERS" + WRITE (UNIT=lunit,FMT="(/,/,T2,A,/)") "SCF PARAMETERS" - WRITE (lunit,"(T3,A,/,T3,A,I5,/,T3,A,ES9.2,/,T3,A,F5.2)")& + WRITE (UNIT=lunit,FMT="(T3,A,/,T3,A,I5,2(/,T3,A,ES9.2),/,T3,A,F5.2)")& "density guess: "//TRIM(density_guess),& "max_scf: ",max_scf,& "eps_scf: ",eps_scf,& + "eps_eigval: ",eps_eigval,& "p_mix: ",p_mix IF (max_diis > 0) THEN - WRITE (lunit,"(T3,A,I5,/,T3,A,ES9.2)")& + WRITE (UNIT=lunit,FMT="(T3,A,I5,/,T3,A,ES9.2)")& "max_diis: ",max_diis,& "eps_diis: ",eps_diis END IF diff --git a/src/radial_solver.F b/src/radial_solver.F index bdf4f5e231..70f098bd98 100644 --- a/src/radial_solver.F +++ b/src/radial_solver.F @@ -96,7 +96,7 @@ CONTAINS ab ( 4, i ) = ab ( 4, i ) * ww END DO - CALL DGBSV ( n, 1, 1, 1, ab, 4, ipiv, s, n, info ) + CALL dgbsv ( n, 1, 1, 1, ab, 4, ipiv, s, n, info ) IF ( info /= 0 ) call stop_program ( "numerov", "DGBSV: info /= 0" ) g ( 1:n ) = s ( 1:n ) diff --git a/src/start_program_run.F b/src/start_program_run.F index 1701727ebc..00dfc073da 100644 --- a/src/start_program_run.F +++ b/src/start_program_run.F @@ -186,7 +186,7 @@ CONTAINS ! *** Local variables *** CHARACTER(LEN=40) :: keyword,test_result - INTEGER :: ipos + INTEGER :: ipos1,ipos2 LOGICAL :: found,print_request ! --------------------------------------------------------------------------- @@ -237,14 +237,16 @@ CONTAINS CALL uppercase(keyword) IF (keyword(1:3) == "NO_") THEN - ipos = 4 + ipos1= 4 print_request = .FALSE. ELSE - ipos = 1 + ipos1= 1 print_request = .TRUE. END IF - SELECT CASE (TRIM(keyword(ipos:))) + ipos2 = LEN_TRIM(keyword) + + SELECT CASE (keyword(ipos1:ipos2)) CASE ("ANGLES") globenv%print%angles = print_request CASE ("ATOMIC_COORDINATES","COORDINATES","COORD") @@ -253,6 +255,8 @@ CONTAINS globenv%print%basic_data_types = print_request CASE ("BASIS_SETS","BASIS_SET","BASIS") globenv%print%basis_set = print_request + CASE ("BLACS_INFORMATION","BLACS_INFO") + globenv%print%blacs_info = print_request CASE ("CARTESIAN_OPERATOR_MATRICES","CARTESIAN_MATRICES") globenv%print%cartesian_om = print_request CASE ("CELL_PARAMETERS","CELL") @@ -408,7 +412,7 @@ CONTAINS "** **",& "** ... make the atoms dance **",& "** **",& - "** Version 3.0 (May 2001) **",& + "** Version 3.0 (June 2001) **",& "** **",& "** Copyright (C) by MPI fuer Festkoerperforschung, Stuttgart (2000) **",& "** **",& diff --git a/src/termination.F b/src/termination.F index 8bca96a507..b3040b7abf 100644 --- a/src/termination.F +++ b/src/termination.F @@ -28,6 +28,7 @@ MODULE termination + USE global_types, ONLY: global_environment_type USE message_passing, ONLY: mp_stop USE output_utilities, ONLY: print_message USE string_utilities, ONLY: compress @@ -223,16 +224,30 @@ CONTAINS !! SOURCE !****************************************************************************** - SUBROUTINE stop_program(routine,error_message) + SUBROUTINE stop_program(routine,error_message,globenv) CHARACTER(LEN=*), INTENT(IN) :: error_message,routine + TYPE(global_environment_type), OPTIONAL, INTENT(IN) :: globenv + +! *** Local variables *** + + LOGICAL :: ionode + ! --------------------------------------------------------------------------- + IF (PRESENT(globenv)) THEN + ionode = globenv%ionode + ELSE + ionode = .TRUE. + END IF + ! *** Print the error message *** - CALL print_message("ERROR in "//TRIM(routine),output_unit,2,2,0) - CALL print_message(error_message,output_unit,1,1,1) + IF (ionode) THEN + CALL print_message("ERROR in "//TRIM(routine),output_unit,2,2,0) + CALL print_message(error_message,output_unit,1,1,1) + END IF CALL mp_stop(0,"stop_program")