Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief various cholesky decomposition related routines
10 : !> \par History
11 : !> 12.2002 Moved routines from cp_cfm_basic_linalg to this new module [Rocco Meli]
12 : ! **************************************************************************************************
13 : MODULE cp_cfm_cholesky
14 : USE cp_cfm_types, ONLY: cp_cfm_type
15 : USE kinds, ONLY: dp
16 : #if defined(__DLAF)
17 : USE cp_cfm_dlaf_api, ONLY: cp_cfm_pzpotrf_dlaf, cp_cfm_pzpotri_dlaf
18 : USE cp_dlaf_utils_api, ONLY: cp_dlaf_initialize, cp_dlaf_create_grid
19 : USE cp_fm_cholesky, ONLY: cholesky_type, dlaf_cholesky_n_min, FM_CHOLESKY_TYPE_DLAF
20 : #endif
21 :
22 : #include "../base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 : PRIVATE
26 :
27 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
28 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_cholesky'
29 :
30 : PUBLIC :: cp_cfm_cholesky_decompose, &
31 : cp_cfm_cholesky_invert
32 :
33 : ! **************************************************************************************************
34 :
35 : CONTAINS
36 :
37 : ! **************************************************************************************************
38 : !> \brief Used to replace a symmetric positive definite matrix M with its Cholesky
39 : !> decomposition U: M = U^T * U, with U upper triangular.
40 : !> \param matrix the matrix to replace with its Cholesky decomposition
41 : !> \param n the number of row (and columns) of the matrix &
42 : !> (defaults to the min(size(matrix)))
43 : !> \param info_out if present, outputs info from (p)zpotrf
44 : !> \par History
45 : !> 05.2002 created [JVdV]
46 : !> 12.2002 updated, added n optional parm [fawzi]
47 : !> 09.2021 removed CPASSERT(info == 0) since there is already check of info [Jan Wilhelm]
48 : !> 12.2024 Added DLA-Future support [Rocco Meli]
49 : !> 03.2025 Moved DLA-Future support to cp_cfm_cholesky_decompose_dlaf [Rocco Meli]
50 : !> \author Joost
51 : ! **************************************************************************************************
52 22356 : SUBROUTINE cp_cfm_cholesky_decompose(matrix, n, info_out)
53 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
54 : INTEGER, INTENT(in), OPTIONAL :: n
55 : INTEGER, INTENT(out), OPTIONAL :: info_out
56 :
57 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_cholesky_decompose'
58 :
59 22356 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
60 : INTEGER :: handle, info, my_n
61 : #if defined(__parallel)
62 : INTEGER, DIMENSION(9) :: desca
63 : #else
64 : INTEGER :: lda
65 : #endif
66 :
67 22356 : CALL timeset(routineN, handle)
68 :
69 : my_n = MIN(matrix%matrix_struct%nrow_global, &
70 22356 : matrix%matrix_struct%ncol_global)
71 22356 : IF (PRESENT(n)) THEN
72 4506 : CPASSERT(n <= my_n)
73 4506 : my_n = n
74 : END IF
75 :
76 22356 : a => matrix%local_data
77 :
78 : #if defined(__parallel)
79 223560 : desca(:) = matrix%matrix_struct%descriptor(:)
80 : #if defined(__DLAF)
81 : IF (cholesky_type == FM_CHOLESKY_TYPE_DLAF .AND. matrix%matrix_struct%nrow_global .GE. dlaf_cholesky_n_min) THEN
82 : ! Initialize DLA-Future on-demand; if already initialized, does nothing
83 : CALL cp_dlaf_initialize()
84 :
85 : ! Create DLAF grid from BLACS context; if already present, does nothing
86 : CALL cp_dlaf_create_grid(matrix%matrix_struct%context%get_handle())
87 :
88 : CALL cp_cfm_pzpotrf_dlaf('U', my_n, a(:, :), 1, 1, desca, info)
89 : ELSE
90 : #endif
91 22356 : CALL pzpotrf('U', my_n, a(1, 1), 1, 1, desca, info)
92 : #if defined(__DLAF)
93 : END IF
94 : #endif
95 :
96 : #else
97 : lda = SIZE(a, 1)
98 : CALL zpotrf('U', my_n, a(1, 1), lda, info)
99 : #endif
100 :
101 22356 : IF (PRESENT(info_out)) THEN
102 4506 : info_out = info
103 : ELSE
104 17850 : IF (info /= 0) &
105 : CALL cp_abort(__LOCATION__, &
106 0 : "Cholesky decompose failed: matrix is not positive definite or ill-conditioned")
107 : END IF
108 :
109 22356 : CALL timestop(handle)
110 :
111 22356 : END SUBROUTINE cp_cfm_cholesky_decompose
112 :
113 : ! **************************************************************************************************
114 : !> \brief Used to replace Cholesky decomposition by the inverse.
115 : !> \param matrix : the matrix to invert (must be an upper triangular matrix),
116 : !> and is the output of Cholesky decomposition
117 : !> \param n : size of the matrix to invert (defaults to the min(size(matrix)))
118 : !> \param info_out : if present, outputs info of (p)zpotri
119 : !> \par History
120 : !> 05.2002 created Lianheng Tong, based on cp_fm_cholesky_invert
121 : !> \author Lianheng Tong
122 : ! **************************************************************************************************
123 4818 : SUBROUTINE cp_cfm_cholesky_invert(matrix, n, info_out)
124 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
125 : INTEGER, INTENT(in), OPTIONAL :: n
126 : INTEGER, INTENT(out), OPTIONAL :: info_out
127 :
128 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_cholesky_invert'
129 4818 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa
130 : INTEGER :: info, handle
131 : INTEGER :: my_n
132 : #if defined(__parallel)
133 : INTEGER, DIMENSION(9) :: desca
134 : #endif
135 :
136 4818 : CALL timeset(routineN, handle)
137 :
138 : my_n = MIN(matrix%matrix_struct%nrow_global, &
139 4818 : matrix%matrix_struct%ncol_global)
140 4818 : IF (PRESENT(n)) THEN
141 0 : CPASSERT(n <= my_n)
142 0 : my_n = n
143 : END IF
144 :
145 4818 : aa => matrix%local_data
146 :
147 : #if defined(__parallel)
148 48180 : desca = matrix%matrix_struct%descriptor
149 : #if defined(__DLAF)
150 : IF (cholesky_type == FM_CHOLESKY_TYPE_DLAF .AND. matrix%matrix_struct%nrow_global .GE. dlaf_cholesky_n_min) THEN
151 : ! Initialize DLA-Future on-demand; if already initialized, does nothing
152 : CALL cp_dlaf_initialize()
153 :
154 : ! Create DLAF grid from BLACS context; if already present, does nothing
155 : CALL cp_dlaf_create_grid(matrix%matrix_struct%context%get_handle())
156 :
157 : CALL cp_cfm_pzpotri_dlaf('U', my_n, aa(:, :), 1, 1, desca, info)
158 : ELSE
159 : #endif
160 4818 : CALL pzpotri('U', my_n, aa(1, 1), 1, 1, desca, info)
161 : #if defined(__DLAF)
162 : END IF
163 : #endif
164 : #else
165 : CALL zpotri('U', my_n, aa(1, 1), SIZE(aa, 1), info)
166 : #endif
167 :
168 4818 : IF (PRESENT(info_out)) THEN
169 0 : info_out = info
170 : ELSE
171 4818 : IF (info /= 0) &
172 : CALL cp_abort(__LOCATION__, &
173 0 : "Cholesky invert failed: the matrix is not positive definite or ill-conditioned.")
174 : END IF
175 :
176 4818 : CALL timestop(handle)
177 :
178 4818 : END SUBROUTINE cp_cfm_cholesky_invert
179 :
180 : END MODULE cp_cfm_cholesky
|
