comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-21 14:35:15 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

- Parallelization of the QS SCF using ScaLAPACK/PBLAS/BLACS (-> blacs.F)

Browse source 
- CRAY and Linux architecture files updated
- Some subroutine names changed to lower case (third try)


svn-origin-rev: 37
This commit is contained in:
Matthias Krack 2001-07-02 14:33:48 +00:00
parent 9bb06770e5
commit 9a894e2d03
28 changed files with 3141 additions and 663 deletions
Download patch file Download diff file Expand all files Collapse all files

25
arch/CRAY-T3E.pdbg
View file

@ -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

25
arch/CRAY-T3E.popt
View file

@ -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

22
arch/CRAY-T3E.sdbg
View file

@ -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

22
arch/CRAY-T3E.sopt
View file

@ -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

15
arch/Linux-i686-pgi.pdbg Normal file
View file

@ -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

15
arch/Linux-i686-pgi.popt Normal file
View file

@ -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

5
arch/Linux-i686-pgi.sdbg
View file

@ -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

5
arch/Linux-i686-pgi.sopt
View file

@ -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

2
src/OBJECTDEFS
View file

@ -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\

1958
src/blacs.F Normal file
View file

File diff suppressed because it is too large Load diff

60
src/collocate_density.F
View file

@ -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

27
src/core_energies.F
View file

@ -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

82
src/core_hamiltonian.F
View file

@ -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

171
src/diis.F
View file

@ -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

17
src/integrate_potential.F
View file

@ -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

12
src/lib/fast.F
View file

@ -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

20
src/lib/mltfftsg.F
View file

@ -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

10
src/library_tests.F
View file

@ -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

4
src/mathlib.F
View file

@ -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 ***

251
src/matrix_types.F
View file

@ -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

57
src/message_passing.F
View file

@ -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
!******************************************************************************

146
src/mo_types.F
View file

@ -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
! *****************************************************************************

250
src/om_utilities.F
View file

@ -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

2
src/print_keys.F
View file

@ -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/))

564
src/qs_scf.F
View file

File diff suppressed because it is too large Load diff

2
src/radial_solver.F
View file

@ -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 )

14
src/start_program_run.F
View file

@ -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) **",&
"** **",&

21
src/termination.F
View file

@ -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")