Skip to content

Commit

Permalink
Two new subroutines which give correct determinant signs even in para…
Browse files Browse the repository at this point in the history 
…llel environment with any number of ranks and threads

1) cp_fm_det computes the determinant of a real square matrix
2) cp_cfm_det computes the determinant of a complex square matrix
Both routines make a copy of the input matrix, thus it is NOT modified in any way and can be safely used after.
Also replaced the calls in two subroutines where determinant sign matters: calc_et_coupling and qs_moment_berry_phase.
  • Loading branch information
@andsinya @mkrack
andsinya authored and mkrack committed Apr 23, 2024
1 parent adfb00c commit b813477
Show file tree
Hide file tree
Showing 4 changed files with 163 additions and 9 deletions.
8 changes: 5 additions & 3 deletions src/et_coupling.F
View file
Expand Up @@ -15,7 +15,8 @@ MODULE et_coupling
USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
dbcsr_deallocate_matrix_set
USE cp_fm_basic_linalg, ONLY: cp_fm_invert,&
cp_fm_transpose
cp_fm_transpose,&
cp_fm_det
USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
fm_pool_get_el_struct
USE cp_fm_struct, ONLY: cp_fm_struct_type
Expand Down Expand Up @@ -156,8 +157,9 @@ SUBROUTINE calc_et_coupling(qs_env)
CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, &
qs_env%et_coupling%et_mo_coeff(i), &
tmp2, 0.0_dp, rest_MO(2))

CALL cp_fm_invert(SMO, inverse_mat, S_det(i))

CALL cp_fm_det(SMO, S_det(i))
CALL cp_fm_invert(SMO, inverse_mat)

CALL cp_fm_get_info(inverse_mat, nrow_local=nrow_local, ncol_local=ncol_local, &
row_indices=row_indices, col_indices=col_indices)
Expand Down
88 changes: 86 additions & 2 deletions src/fm/cp_cfm_basic_linalg.F
View file
Expand Up @@ -16,7 +16,11 @@
MODULE cp_cfm_basic_linalg
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_cfm_types, ONLY: cp_cfm_get_info,&
cp_cfm_type
cp_cfm_type,&
cp_cfm_to_cfm,&
cp_cfm_create,&
cp_cfm_release
USE cp_log_handling, ONLY: cp_to_string
USE cp_fm_struct, ONLY: cp_fm_struct_equivalent
USE cp_fm_types, ONLY: cp_fm_type
USE kahan_sum, ONLY: accurate_dot_product
Expand Down Expand Up @@ -49,7 +53,8 @@ MODULE cp_cfm_basic_linalg
cp_cfm_triangular_multiply, &
cp_cfm_rot_rows, &
cp_cfm_rot_cols, &
cp_cfm_cholesky_invert
cp_cfm_cholesky_invert, &
cp_cfm_det ! determinant of a complex matrix with correct sign
REAL(kind=dp), EXTERNAL :: zlange, pzlange
Expand All @@ -61,6 +66,85 @@ MODULE cp_cfm_basic_linalg
CONTAINS
! **************************************************************************************************
!> \brief Computes the determinant (with a correct sign even in parallel environment!) of a complex square matrix
!> \author A. Sinyavskiy (andrey.sinyavskiy@chem.uzh.ch)
! **************************************************************************************************
SUBROUTINE cp_cfm_det(matrix_a, det_a)
TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
COMPLEX(KIND=dp), INTENT(OUT) :: det_a
COMPLEX(KIND=dp) :: determinant
TYPE(cp_cfm_type) :: matrix_lu
COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a
INTEGER :: n, i, info, P
INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
COMPLEX(KIND=dp), DIMENSION(:), POINTER :: diag
#if defined(__SCALAPACK)
INTEGER, EXTERNAL :: indxl2g
INTEGER :: myprow, nprow, npcol, nrow_local, nrow_block, irow_local, &
mypcol, ncol_local, ncol_block, icol_local, j
INTEGER, DIMENSION(9) :: desca
#endif
CALL cp_cfm_create(matrix=matrix_lu, &
matrix_struct=matrix_a%matrix_struct, &
name="A_lu"//TRIM(ADJUSTL(cp_to_string(1)))//"MATRIX")
CALL cp_cfm_to_cfm(matrix_a, matrix_lu)
a => matrix_lu%local_data
n = matrix_lu%matrix_struct%nrow_global
ALLOCATE (ipivot(n))
ipivot(:) = 0
P = 0
ALLOCATE (diag(n))
diag(:) = 0.0_dp
#if defined(__SCALAPACK)
! Use LU decomposition
desca(:) = matrix_lu%matrix_struct%descriptor(:)
CALL pzgetrf(n, n, a(1, 1), 1, 1, desca, ipivot, info)
myprow = matrix_lu%matrix_struct%context%mepos(1)
mypcol = matrix_lu%matrix_struct%context%mepos(2)
nprow = matrix_lu%matrix_struct%context%num_pe(1)
npcol = matrix_lu%matrix_struct%context%num_pe(2)
nrow_local = matrix_lu%matrix_struct%nrow_locals(myprow)
nrow_block = matrix_lu%matrix_struct%nrow_block
ncol_block = matrix_lu%matrix_struct%ncol_block
ncol_local = matrix_lu%matrix_struct%ncol_locals(mypcol)
DO irow_local = 1, nrow_local
i = indxl2g(irow_local, nrow_block, myprow, matrix_lu%matrix_struct%first_p_pos(1), nprow)
DO icol_local = 1, ncol_local
j = indxl2g(icol_local, ncol_block, mypcol, matrix_lu%matrix_struct%first_p_pos(2), npcol)
IF (i == j) diag(i) = matrix_lu%local_data(irow_local, icol_local)
END DO
END DO
CALL matrix_lu%matrix_struct%para_env%sum(diag)
determinant = PRODUCT(diag)
DO irow_local = 1, nrow_local
i = indxl2g(irow_local, nrow_block, myprow, matrix_lu%matrix_struct%first_p_pos(1), nprow)
IF (ipivot(irow_local) /= i) P = P + 1
END DO
CALL matrix_lu%matrix_struct%para_env%sum(P)
! very important fix
P = P/npcol
#else
CALL zgetrf(n, n, a(1, 1), n, ipivot, info)
DO i = 1, n
diag(i) = matrix_lu%local_data(i, i)
END DO
determinant = PRODUCT(diag)
DO i = 1, n
IF (ipivot(i) /= i) P = P + 1
END DO
#endif
DEALLOCATE (ipivot)
DEALLOCATE (diag)
CALL cp_cfm_release(matrix_lu)
det_a = determinant*(-2*MOD(P, 2) + 1.0_dp)
END SUBROUTINE cp_cfm_det
! **************************************************************************************************
!> \brief Computes the element-wise (Schur) product of two matrices: C = A \circ B .
!> \param matrix_a the first input matrix
Expand Down
68 changes: 67 additions & 1 deletion src/fm/cp_fm_basic_linalg.F
View file
Expand Up @@ -65,7 +65,8 @@ MODULE cp_fm_basic_linalg
cp_fm_rot_rows, & ! rotates two rows
cp_fm_rot_cols, & ! rotates two columns
cp_fm_cholesky_restore, & ! apply Cholesky decomposition
cp_fm_Gram_Schmidt_orthonorm ! Gram-Schmidt orthonormalization of columns of a full matrix
cp_fm_Gram_Schmidt_orthonorm, & ! Gram-Schmidt orthonormalization of columns of a full matrix, &
cp_fm_det ! determinant of a real matrix with correct sign

REAL(KIND=dp), EXTERNAL :: dlange, pdlange, pdlatra
REAL(KIND=sp), EXTERNAL :: slange, pslange, pslatra
Expand All @@ -88,6 +89,71 @@ MODULE cp_fm_basic_linalg
END INTERFACE cp_fm_contracted_trace
CONTAINS

! **************************************************************************************************
!> \brief Computes the determinant (with a correct sign even in parallel environment!) of a real square matrix
!> \author A. Sinyavskiy (andrey.sinyavskiy@chem.uzh.ch)
! **************************************************************************************************
SUBROUTINE cp_fm_det(matrix_a, det_a)

TYPE(cp_fm_type), INTENT(IN) :: matrix_a
REAL(KIND=dp), INTENT(OUT) :: det_a
REAL(KIND=dp) :: determinant
TYPE(cp_fm_type) :: matrix_lu
REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
INTEGER :: n, i, info, P
INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
REAL(KIND=dp), DIMENSION(:), POINTER :: diag

#if defined(__SCALAPACK)
INTEGER, EXTERNAL :: indxl2g
INTEGER :: myprow, nprow, npcol, nrow_local, nrow_block, irow_local
INTEGER, DIMENSION(9) :: desca
#endif

CALL cp_fm_create(matrix=matrix_lu, &
matrix_struct=matrix_a%matrix_struct, &
name="A_lu"//TRIM(ADJUSTL(cp_to_string(1)))//"MATRIX")
CALL cp_fm_to_fm(matrix_a, matrix_lu)

a => matrix_lu%local_data
n = matrix_lu%matrix_struct%nrow_global
ALLOCATE (ipivot(n))
ipivot(:) = 0
P = 0
ALLOCATE (diag(n))
diag(:) = 0.0_dp
#if defined(__SCALAPACK)
! Use LU decomposition
desca(:) = matrix_lu%matrix_struct%descriptor(:)
CALL pdgetrf(n, n, a, 1, 1, desca, ipivot, info)
CALL cp_fm_get_diag(matrix_lu, diag)
determinant = PRODUCT(diag)
myprow = matrix_lu%matrix_struct%context%mepos(1)
nprow = matrix_lu%matrix_struct%context%num_pe(1)
npcol = matrix_lu%matrix_struct%context%num_pe(2)
nrow_local = matrix_lu%matrix_struct%nrow_locals(myprow)
nrow_block = matrix_lu%matrix_struct%nrow_block
DO irow_local = 1, nrow_local
i = indxl2g(irow_local, nrow_block, myprow, matrix_lu%matrix_struct%first_p_pos(1), nprow)
IF (ipivot(irow_local) /= i) P = P + 1
END DO
CALL matrix_lu%matrix_struct%para_env%sum(P)
! very important fix
P = P/npcol
#else
CALL dgetrf(n, n, a, n, ipivot, info)
CALL cp_fm_get_diag(matrix_lu, diag)
determinant = PRODUCT(diag)
DO i = 1, n
IF (ipivot(i) /= i) P = P + 1
END DO
#endif
DEALLOCATE (ipivot)
DEALLOCATE (diag)
CALL cp_fm_release(matrix_lu)
det_a = determinant*(-2*MOD(P, 2) + 1.0_dp)
END SUBROUTINE cp_fm_det

! **************************************************************************************************
!> \brief calc A <- alpha*A + beta*B
!> optimized for alpha == 1.0 (just add beta*B) and beta == 0.0 (just
Expand Down
8 changes: 5 additions & 3 deletions src/qs_moments.F
View file
Expand Up @@ -29,7 +29,7 @@ MODULE qs_moments
USE cell_types, ONLY: cell_type,&
pbc
USE commutator_rpnl, ONLY: build_com_mom_nl
USE cp_cfm_basic_linalg, ONLY: cp_cfm_lu_decompose
USE cp_cfm_basic_linalg, ONLY: cp_cfm_det
USE cp_cfm_types, ONLY: cp_cfm_create,&
cp_cfm_get_info,&
cp_cfm_release,&
Expand Down Expand Up @@ -1971,7 +1971,8 @@ SUBROUTINE qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_poin
CMPLX(op_fm_set(1, ispin)%local_data(:, idim), &
-op_fm_set(2, ispin)%local_data(:, idim), dp)
END DO
CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
CALL cp_cfm_det(eigrmat(ispin), zdeta)
zdet = zdet*zdeta
IF (dft_control%nspins == 1) THEN
zdet = zdet*zdeta
Expand Down Expand Up @@ -2002,7 +2003,8 @@ SUBROUTINE qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_poin
CMPLX(op_fm_set(1, ispin)%local_data(:, idim), &
-op_fm_set(2, ispin)%local_data(:, idim), dp)
END DO
CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
CALL cp_cfm_det(eigrmat(ispin), zdeta)
zdet = zdet*zdeta
IF (dft_control%nspins == 1) THEN
zdet = zdet*zdeta
Expand Down

0 comments on commit b813477

Please sign in to comment.