From c08dc9e0b62b5cc4b7a3e2d1c0cfcbdce7cf24e8 Mon Sep 17 00:00:00 2001 From: Hans Pabst Date: Thu, 14 Aug 2025 10:48:17 +0200 Subject: [PATCH] Fixes and out-of-place trsm/trmm - Removed cp_fm_cholesky_restore from cp_fm_basic_linalg - Thereby reduced duplicated code (cp_fm_cholesky_restore). - Prefer cp_fm_cholesky to provide this routine. - Fixed cp_fm_cholesky_restore for pos=right. - Fixed cp_fm_triangular_multiply for side=R. - Removed unused variables. - Adjusted some INTENTs. --- src/base/base_uses.f90 | 7 ++ src/cp_dbcsr_cholesky.F | 4 +- src/fm/cp_fm_basic_linalg.F | 110 ++---------------------- src/fm/cp_fm_cholesky.F | 161 +++++++++++++++++------------------- 4 files changed, 93 insertions(+), 189 deletions(-) diff --git a/src/base/base_uses.f90 b/src/base/base_uses.f90 index abbbdb1b08..9378b131fc 100644 --- a/src/base/base_uses.f90 +++ b/src/base/base_uses.f90 @@ -76,3 +76,10 @@ #else #define OMP_DEFAULT_NONE_WITH_OOP NONE #endif + +#if defined (__MKL) +! Preprocessing is enabled by default, and below header is not language specific +! Defines __INTEL_MKL__ (2025), __INTEL_MKL_MINOR__ (0), __INTEL_MKL_UPDATE__ (2), +! and __INTEL_MKL_PATCH__ (0) as well INTEL_MKL_VERSION, e.g., 20250200. +#include +#endif diff --git a/src/cp_dbcsr_cholesky.F b/src/cp_dbcsr_cholesky.F index 4cdabaf911..a9fcc392a2 100644 --- a/src/cp_dbcsr_cholesky.F +++ b/src/cp_dbcsr_cholesky.F @@ -20,10 +20,10 @@ MODULE cp_dbcsr_cholesky dbcsr_type USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& copy_fm_to_dbcsr - USE cp_fm_basic_linalg, ONLY: cp_fm_cholesky_restore,& - cp_fm_potrf,& + USE cp_fm_basic_linalg, ONLY: cp_fm_potrf,& cp_fm_potri,& cp_fm_uplo_to_full + USE cp_fm_cholesky, ONLY: cp_fm_cholesky_restore USE cp_fm_struct, ONLY: cp_fm_struct_create,& cp_fm_struct_release,& cp_fm_struct_type diff --git a/src/fm/cp_fm_basic_linalg.F b/src/fm/cp_fm_basic_linalg.F index f563a59782..1f2fd19779 100644 --- a/src/fm/cp_fm_basic_linalg.F +++ b/src/fm/cp_fm_basic_linalg.F @@ -64,7 +64,6 @@ MODULE cp_fm_basic_linalg cp_fm_potri, & ! Invert triangular matrix 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_det, & ! determinant of a real matrix with correct sign cp_fm_matvec ! matrix-vector multiplication (vector replicated) @@ -716,7 +715,7 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_schur_product' INTEGER :: handle, icol_local, irow_local, mypcol, & - myprow, ncol_local, npcol, nprow, & + myprow, ncol_local, & nrow_local REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, c TYPE(cp_blacs_env_type), POINTER :: context @@ -726,8 +725,6 @@ CONTAINS context => matrix_a%matrix_struct%context myprow = context%mepos(1) mypcol = context%mepos(2) - nprow = context%num_pe(1) - npcol = context%num_pe(2) a => matrix_a%local_data b => matrix_b%local_data @@ -767,7 +764,7 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_trace_a0b0t0' INTEGER :: handle, mypcol, myprow, ncol_local, & - npcol, nprow, nrow_local + nrow_local REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp, b_sp TYPE(cp_blacs_env_type), POINTER :: context @@ -778,8 +775,6 @@ CONTAINS context => matrix_a%matrix_struct%context myprow = context%mepos(1) mypcol = context%mepos(2) - nprow = context%num_pe(1) - npcol = context%num_pe(2) group = matrix_a%matrix_struct%para_env @@ -1117,7 +1112,7 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_triangular_multiply' CHARACTER :: side_char, transa, unit_diag, uplo - INTEGER :: handle, m, n + INTEGER :: handle, mdim, m, n LOGICAL :: invert REAL(KIND=dp) :: al @@ -1130,6 +1125,7 @@ CONTAINS al = 1.0_dp CALL cp_fm_get_info(matrix_b, nrow_global=m, ncol_global=n) IF (PRESENT(side)) side_char = side + mdim = MERGE(1, 2, 'L' .EQ. side_char) IF (PRESENT(invert_tr)) invert = invert_tr IF (PRESENT(uplo_tr)) uplo = uplo_tr IF (PRESENT(unit_diag_tr)) THEN @@ -1161,7 +1157,7 @@ CONTAINS #else CALL dtrsm(side_char, uplo, transa, unit_diag, m, n, al, & triangular_matrix%local_data(1, 1), & - SIZE(triangular_matrix%local_data, 1), & + SIZE(triangular_matrix%local_data, mdim), & matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1)) #endif @@ -1176,7 +1172,7 @@ CONTAINS #else CALL dtrmm(side_char, uplo, transa, unit_diag, m, n, al, & triangular_matrix%local_data(1, 1), & - SIZE(triangular_matrix%local_data, 1), & + SIZE(triangular_matrix%local_data, mdim), & matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1)) #endif @@ -1285,7 +1281,7 @@ CONTAINS CHARACTER :: myuplo INTEGER :: handle, icol_global, irow_global, & mypcol, myprow, ncol_global, & - npcol, nprow, nrow_global + nrow_global REAL(KIND=dp), DIMENSION(:, :), POINTER :: a REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp TYPE(cp_blacs_env_type), POINTER :: context @@ -1315,8 +1311,6 @@ CONTAINS context => matrix%matrix_struct%context myprow = context%mepos(1) mypcol = context%mepos(2) - nprow = context%num_pe(1) - npcol = context%num_pe(2) #if defined(__parallel) @@ -2711,96 +2705,6 @@ CONTAINS CPASSERT(info == 0) END SUBROUTINE cp_fm_potri -! ************************************************************************************************** -!> \brief ... -!> \param fm_matrix ... -!> \param neig ... -!> \param fm_matrixb ... -!> \param fm_matrixout ... -!> \param op ... -!> \param pos ... -!> \param transa ... -! ************************************************************************************************** - SUBROUTINE cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa) - TYPE(cp_fm_type) :: fm_matrix - TYPE(cp_fm_type) :: fm_matrixb - TYPE(cp_fm_type) :: fm_matrixout - INTEGER, INTENT(IN) :: neig - CHARACTER(LEN=*), INTENT(IN) :: op - CHARACTER(LEN=*), INTENT(IN) :: pos - CHARACTER(LEN=*), INTENT(IN) :: transa - - REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, outm - REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp, b_sp, outm_sp - INTEGER :: n, itype - REAL(KIND=dp) :: alpha -#if defined(__parallel) - INTEGER :: i - INTEGER, DIMENSION(9) :: desca, descb, descout -#endif - - ! notice b is the cholesky guy - a => fm_matrix%local_data - b => fm_matrixb%local_data - outm => fm_matrixout%local_data - a_sp => fm_matrix%local_data_sp - b_sp => fm_matrixb%local_data_sp - outm_sp => fm_matrixout%local_data_sp - - n = fm_matrix%matrix_struct%nrow_global - itype = 1 - -#if defined(__parallel) - desca(:) = fm_matrix%matrix_struct%descriptor(:) - descb(:) = fm_matrixb%matrix_struct%descriptor(:) - descout(:) = fm_matrixout%matrix_struct%descriptor(:) - alpha = 1.0_dp - DO i = 1, neig - IF (fm_matrix%use_sp) THEN - CALL pscopy(n, a_sp(1, 1), 1, i, desca, 1, outm_sp(1, 1), 1, i, descout, 1) - ELSE - CALL pdcopy(n, a(1, 1), 1, i, desca, 1, outm(1, 1), 1, i, descout, 1) - END IF - END DO - IF (op .EQ. "SOLVE") THEN - IF (fm_matrix%use_sp) THEN - CALL pstrsm(pos, 'U', transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), 1, 1, descb, & - outm_sp(1, 1), 1, 1, descout) - ELSE - CALL pdtrsm(pos, 'U', transa, 'N', n, neig, alpha, b(1, 1), 1, 1, descb, outm(1, 1), 1, 1, descout) - END IF - ELSE - IF (fm_matrix%use_sp) THEN - CALL pstrmm(pos, 'U', transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), 1, 1, descb, & - outm_sp(1, 1), 1, 1, descout) - ELSE - CALL pdtrmm(pos, 'U', transa, 'N', n, neig, alpha, b(1, 1), 1, 1, descb, outm(1, 1), 1, 1, descout) - END IF - END IF -#else - alpha = 1.0_dp - IF (fm_matrix%use_sp) THEN - CALL scopy(neig*n, a_sp(1, 1), 1, outm_sp(1, 1), 1) - ELSE - CALL dcopy(neig*n, a(1, 1), 1, outm(1, 1), 1) - END IF - IF (op .EQ. "SOLVE") THEN - IF (fm_matrix%use_sp) THEN - CALL strsm(pos, 'U', transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), SIZE(b_sp, 1), outm_sp(1, 1), n) - ELSE - CALL dtrsm(pos, 'U', transa, 'N', n, neig, alpha, b(1, 1), SIZE(b, 1), outm(1, 1), n) - END IF - ELSE - IF (fm_matrix%use_sp) THEN - CALL strmm(pos, 'U', transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), n, outm_sp(1, 1), n) - ELSE - CALL dtrmm(pos, 'U', transa, 'N', n, neig, alpha, b(1, 1), n, outm(1, 1), n) - END IF - END IF -#endif - - END SUBROUTINE cp_fm_cholesky_restore - ! ************************************************************************************************** !> \brief Calculates !> yv = alpha*amat*xv + beta*yv diff --git a/src/fm/cp_fm_cholesky.F b/src/fm/cp_fm_cholesky.F index 357c754d08..67d3a8d6c5 100644 --- a/src/fm/cp_fm_cholesky.F +++ b/src/fm/cp_fm_cholesky.F @@ -12,7 +12,6 @@ !> \author Fawzi Mohamed ! ************************************************************************************************** MODULE cp_fm_cholesky - USE cp_blacs_env, ONLY: cp_blacs_env_type USE cp_dlaf_utils_api, ONLY: cp_dlaf_create_grid,& cp_dlaf_initialize USE cp_fm_dlaf_api, ONLY: cp_pdpotrf_dlaf,& @@ -140,13 +139,13 @@ CONTAINS !> \author Joost VandeVondele ! ************************************************************************************************** SUBROUTINE cp_fm_cholesky_invert(matrix, n, info_out) - TYPE(cp_fm_type), INTENT(IN) :: matrix + TYPE(cp_fm_type), INTENT(IN) :: matrix INTEGER, INTENT(in), OPTIONAL :: n INTEGER, INTENT(OUT), OPTIONAL :: info_out CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_invert' - REAL(KIND=dp), DIMENSION(:, :), POINTER :: a - REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp + REAL(KIND=dp), DIMENSION(:, :), POINTER :: a + REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp INTEGER :: info, handle INTEGER :: my_n #if defined(__parallel) @@ -228,11 +227,11 @@ CONTAINS !> \author Joost VandeVondele ! ************************************************************************************************** SUBROUTINE cp_fm_cholesky_reduce(matrix, matrixb, itype) - TYPE(cp_fm_type), INTENT(IN) :: matrix, matrixb - INTEGER, OPTIONAL :: itype + TYPE(cp_fm_type), INTENT(IN) :: matrix, matrixb + INTEGER, OPTIONAL :: itype CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_reduce' - REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b + REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b INTEGER :: info, handle INTEGER :: n, my_itype #if defined(__parallel) @@ -273,130 +272,124 @@ CONTAINS END SUBROUTINE cp_fm_cholesky_reduce -! -! op can be "SOLVE" (out = U^-1 * in ) or "MULTIPLY" (out = U * in ) -! pos can be "LEFT" or "RIGHT" (U at the left or at the right) -! -! DEPRECATED, see cp_fm_basic_linalg:cp_fm_triangular_multiply -! ! ************************************************************************************************** -!> \brief ... -!> \param matrix ... +!> \brief apply Cholesky decomposition +!> op can be "SOLVE" (out = U^-1 * in) or "MULTIPLY" (out = U * in) +!> pos can be "LEFT" or "RIGHT" (U at the left or at the right) +!> \param fm_matrix ... !> \param neig ... -!> \param matrixb ... -!> \param matrixout ... +!> \param fm_matrixb ... +!> \param fm_matrixout ... !> \param op ... !> \param pos ... !> \param transa ... ! ************************************************************************************************** - SUBROUTINE cp_fm_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa) - TYPE(cp_fm_type), INTENT(IN) :: matrix, matrixb, matrixout - INTEGER, INTENT(IN) :: neig - CHARACTER(LEN=*), INTENT(IN) :: op - CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: pos - CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: transa + SUBROUTINE cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa) + TYPE(cp_fm_type), INTENT(IN) :: fm_matrix, fm_matrixb, fm_matrixout + INTEGER, INTENT(IN) :: neig + CHARACTER(LEN=*), INTENT(IN) :: op + CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: pos + CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: transa CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_restore' - REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, out - REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp, b_sp, out_sp - INTEGER :: itype, handle - INTEGER :: n - REAL(KIND=dp) :: alpha - INTEGER :: myprow, mypcol - TYPE(cp_blacs_env_type), POINTER :: context + REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, outm + REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp, b_sp, outm_sp + REAL(KIND=dp) :: alpha + INTEGER :: handle, mdim, n CHARACTER :: chol_pos, chol_transa #if defined(__parallel) INTEGER :: i INTEGER, DIMENSION(9) :: desca, descb, descout +#elif defined(__MKL) && (2025 <= __INTEL_MKL__) + REAL(KIND=dp) :: beta #endif CALL timeset(routineN, handle) - context => matrix%matrix_struct%context - myprow = context%mepos(1) - mypcol = context%mepos(2) - n = matrix%matrix_struct%nrow_global - itype = 1 - IF (op /= "SOLVE" .AND. op /= "MULTIPLY") & - CPABORT("wrong argument op") + ! wrong argument op + CPASSERT((op == "SOLVE") .OR. (op == "MULTIPLY")) + ! not the same precision + CPASSERT(fm_matrix%use_sp .EQV. fm_matrixb%use_sp) + CPASSERT(fm_matrix%use_sp .EQV. fm_matrixout%use_sp) - IF (PRESENT(pos)) THEN - SELECT CASE (pos) - CASE ("LEFT") - chol_pos = 'L' - CASE ("RIGHT") - chol_pos = 'R' - CASE DEFAULT - CPABORT("wrong argument pos") - END SELECT - ELSE - chol_pos = 'L' - END IF + chol_pos = 'L' + IF (PRESENT(pos)) chol_pos = pos(1:1) + CPASSERT((chol_pos == 'L') .OR. (chol_pos == 'R')) chol_transa = 'N' IF (PRESENT(transa)) chol_transa = transa - IF ((matrix%use_sp .NEQV. matrixb%use_sp) .OR. (matrix%use_sp .NEQV. matrixout%use_sp)) & - CPABORT("not the same precision") - ! notice b is the cholesky guy - a => matrix%local_data - b => matrixb%local_data - out => matrixout%local_data - a_sp => matrix%local_data_sp - b_sp => matrixb%local_data_sp - out_sp => matrixout%local_data_sp + a => fm_matrix%local_data + b => fm_matrixb%local_data + outm => fm_matrixout%local_data + a_sp => fm_matrix%local_data_sp + b_sp => fm_matrixb%local_data_sp + outm_sp => fm_matrixout%local_data_sp + n = fm_matrix%matrix_struct%nrow_global + mdim = MERGE(1, 2, 'L' .EQ. chol_pos) + alpha = 1.0_dp #if defined(__parallel) - - desca(:) = matrix%matrix_struct%descriptor(:) - descb(:) = matrixb%matrix_struct%descriptor(:) - descout(:) = matrixout%matrix_struct%descriptor(:) - alpha = 1.0_dp + desca(:) = fm_matrix%matrix_struct%descriptor(:) + descb(:) = fm_matrixb%matrix_struct%descriptor(:) + descout(:) = fm_matrixout%matrix_struct%descriptor(:) DO i = 1, neig - IF (matrix%use_sp) THEN - CALL pscopy(n, a_sp(1, 1), 1, i, desca, 1, out_sp(1, 1), 1, i, descout, 1) + IF (fm_matrix%use_sp) THEN + CALL pscopy(n, a_sp(1, 1), 1, i, desca, 1, outm_sp(1, 1), 1, i, descout, 1) ELSE - CALL pdcopy(n, a(1, 1), 1, i, desca, 1, out(1, 1), 1, i, descout, 1) + CALL pdcopy(n, a(1, 1), 1, i, desca, 1, outm(1, 1), 1, i, descout, 1) END IF END DO IF (op .EQ. "SOLVE") THEN - IF (matrix%use_sp) THEN + IF (fm_matrix%use_sp) THEN CALL pstrsm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), 1, 1, descb, & - out_sp(1, 1), 1, 1, descout) + outm_sp(1, 1), 1, 1, descout) ELSE - CALL pdtrsm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), 1, 1, descb, out(1, 1), 1, 1, descout) + CALL pdtrsm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), 1, 1, descb, outm(1, 1), 1, 1, descout) END IF ELSE - IF (matrix%use_sp) THEN + IF (fm_matrix%use_sp) THEN CALL pstrmm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), 1, 1, descb, & - out_sp(1, 1), 1, 1, descout) + outm_sp(1, 1), 1, 1, descout) ELSE - CALL pdtrmm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), 1, 1, descb, out(1, 1), 1, 1, descout) + CALL pdtrmm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), 1, 1, descb, outm(1, 1), 1, 1, descout) + END IF + END IF +#elif defined(__MKL) && (2025 <= __INTEL_MKL__) + beta = 0.0_dp + IF (op .EQ. "SOLVE") THEN + IF (fm_matrix%use_sp) THEN + CALL strsm_oop(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), SIZE(b_sp, mdim), a_sp(1, 1), n, REAL(beta, sp), outm_sp(1, 1), n) + ELSE + CALL dtrsm_oop(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), SIZE(b, mdim), a(1, 1), n, beta, outm(1, 1), n) + END IF + ELSE + IF (fm_matrix%use_sp) THEN + CALL strmm_oop(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), SIZE(b_sp, mdim), a_sp(1, 1), n, REAL(beta, sp), outm_sp(1, 1), n) + ELSE + CALL dtrmm_oop(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), SIZE(b, mdim), a(1, 1), n, beta, outm(1, 1), n) END IF END IF #else - - alpha = 1.0_dp - IF (matrix%use_sp) THEN - CALL scopy(neig*n, a_sp(1, 1), 1, out_sp(1, 1), 1) + IF (fm_matrix%use_sp) THEN + CALL scopy(neig*n, a_sp(1, 1), 1, outm_sp(1, 1), 1) ELSE - CALL dcopy(neig*n, a(1, 1), 1, out(1, 1), 1) + CALL dcopy(neig*n, a(1, 1), 1, outm(1, 1), 1) END IF IF (op .EQ. "SOLVE") THEN - IF (matrix%use_sp) THEN - CALL strsm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), SIZE(b_sp, 1), out_sp(1, 1), n) + IF (fm_matrix%use_sp) THEN + CALL strsm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), SIZE(b_sp, mdim), outm_sp(1, 1), n) ELSE - CALL dtrsm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), SIZE(b, 1), out(1, 1), n) + CALL dtrsm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), SIZE(b, mdim), outm(1, 1), n) END IF ELSE - IF (matrix%use_sp) THEN - CALL strmm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), n, out_sp(1, 1), n) + IF (fm_matrix%use_sp) THEN + CALL strmm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), SIZE(b_sp, mdim), outm_sp(1, 1), n) ELSE - CALL dtrmm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), n, out(1, 1), n) + CALL dtrmm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), SIZE(b, mdim), outm(1, 1), n) END IF END IF - #endif CALL timestop(handle)