(git:ccc2433)
cp_dbcsr_cholesky.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Interface to (sca)lapack for the Cholesky based procedures
10 !> \author VW
11 !> \date 2009-09-08
12 !> \version 0.8
13 !>
14 !> <b>Modification history:</b>
15 !> - Created 2009-09-08
16 ! **************************************************************************************************
17 MODULE cp_dbcsr_cholesky
18  USE cp_blacs_env, ONLY: cp_blacs_env_type
19  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
20  copy_fm_to_dbcsr
21  USE cp_fm_basic_linalg, ONLY: cp_fm_cholesky_restore,&
22  cp_fm_potrf,&
23  cp_fm_potri,&
24  cp_fm_upper_to_full
25  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
26  cp_fm_struct_release,&
27  cp_fm_struct_type
28  USE cp_fm_types, ONLY: cp_fm_create,&
29  cp_fm_release,&
30  cp_fm_type
31  USE dbcsr_api, ONLY: dbcsr_get_info,&
32  dbcsr_type
33  USE message_passing, ONLY: mp_para_env_type
34 #include "base/base_uses.f90"
35 
36  IMPLICIT NONE
37 
38  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_dbcsr_cholesky'
39 
40  PUBLIC :: cp_dbcsr_cholesky_decompose, cp_dbcsr_cholesky_invert, &
41  cp_dbcsr_cholesky_restore
42 
43  PRIVATE
44 
45 CONTAINS
46 
47 ! **************************************************************************************************
48 !> \brief used to replace a symmetric positive def. matrix M with its cholesky
49 !> decomposition U: M = U^T * U, with U upper triangular
50 !> \param matrix the matrix to replace with its cholesky decomposition
51 !> \param n the number of row (and columns) of the matrix &
52 !> (defaults to the min(size(matrix)))
53 !> \param para_env ...
54 !> \param blacs_env ...
55 !> \par History
56 !> 05.2002 created [JVdV]
57 !> 12.2002 updated, added n optional parm [fawzi]
58 !> \author Joost
59 ! **************************************************************************************************
60  SUBROUTINE cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
61  TYPE(dbcsr_type) :: matrix
62  INTEGER, INTENT(in), OPTIONAL :: n
63  TYPE(mp_para_env_type), POINTER :: para_env
64  TYPE(cp_blacs_env_type), POINTER :: blacs_env
65 
66  CHARACTER(len=*), PARAMETER :: routinen = 'cp_dbcsr_cholesky_decompose'
67 
68  INTEGER :: handle, my_n, nfullcols_total, &
69  nfullrows_total
70  TYPE(cp_fm_struct_type), POINTER :: fm_struct
71  TYPE(cp_fm_type) :: fm_matrix
72 
73  CALL timeset(routinen, handle)
74 
75  NULLIFY (fm_struct)
76  CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
77 
78  CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
79  ncol_global=nfullcols_total, para_env=para_env)
80  CALL cp_fm_create(fm_matrix, fm_struct, name="fm_matrix")
81  CALL cp_fm_struct_release(fm_struct)
82 
83  CALL copy_dbcsr_to_fm(matrix, fm_matrix)
84 
85  my_n = min(fm_matrix%matrix_struct%nrow_global, &
86  fm_matrix%matrix_struct%ncol_global)
87  IF (PRESENT(n)) THEN
88  cpassert(n <= my_n)
89  my_n = n
90  END IF
91 
92  CALL cp_fm_potrf(fm_matrix, my_n)
93 
94  CALL copy_fm_to_dbcsr(fm_matrix, matrix)
95 
96  CALL cp_fm_release(fm_matrix)
97 
98  CALL timestop(handle)
99 
100  END SUBROUTINE cp_dbcsr_cholesky_decompose
101 
102 ! **************************************************************************************************
103 !> \brief used to replace the cholesky decomposition by the inverse
104 !> \param matrix the matrix to invert (must be an upper triangular matrix)
105 !> \param n size of the matrix to invert (defaults to the min(size(matrix)))
106 !> \param para_env ...
107 !> \param blacs_env ...
108 !> \param upper_to_full ...
109 !> \par History
110 !> 05.2002 created [JVdV]
111 !> \author Joost VandeVondele
112 ! **************************************************************************************************
113  SUBROUTINE cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, upper_to_full)
114  TYPE(dbcsr_type) :: matrix
115  INTEGER, INTENT(in), OPTIONAL :: n
116  TYPE(mp_para_env_type), POINTER :: para_env
117  TYPE(cp_blacs_env_type), POINTER :: blacs_env
118  LOGICAL, INTENT(IN) :: upper_to_full
119 
120  CHARACTER(len=*), PARAMETER :: routinen = 'cp_dbcsr_cholesky_invert'
121 
122  INTEGER :: handle, my_n, nfullcols_total, &
123  nfullrows_total
124  TYPE(cp_fm_struct_type), POINTER :: fm_struct
125  TYPE(cp_fm_type) :: fm_matrix, fm_matrix_tmp
126 
127  CALL timeset(routinen, handle)
128 
129  NULLIFY (fm_struct)
130  CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
131 
132  CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
133  ncol_global=nfullrows_total, para_env=para_env)
134  CALL cp_fm_create(fm_matrix, fm_struct, name="fm_matrix")
135  CALL cp_fm_struct_release(fm_struct)
136 
137  CALL copy_dbcsr_to_fm(matrix, fm_matrix)
138 
139  my_n = min(fm_matrix%matrix_struct%nrow_global, &
140  fm_matrix%matrix_struct%ncol_global)
141  IF (PRESENT(n)) THEN
142  cpassert(n <= my_n)
143  my_n = n
144  END IF
145 
146  CALL cp_fm_potri(fm_matrix, my_n)
147 
148  IF (upper_to_full) THEN
149  CALL cp_fm_create(fm_matrix_tmp, fm_matrix%matrix_struct, name="fm_matrix_tmp")
150  CALL cp_fm_upper_to_full(fm_matrix, fm_matrix_tmp)
151  CALL cp_fm_release(fm_matrix_tmp)
152  END IF
153 
154  CALL copy_fm_to_dbcsr(fm_matrix, matrix)
155 
156  CALL cp_fm_release(fm_matrix)
157 
158  CALL timestop(handle)
159 
160  END SUBROUTINE cp_dbcsr_cholesky_invert
161 
162 ! **************************************************************************************************
163 !> \brief ...
164 !> \param matrix ...
165 !> \param neig ...
166 !> \param matrixb ...
167 !> \param matrixout ...
168 !> \param op ...
169 !> \param pos ...
170 !> \param transa ...
171 !> \param para_env ...
172 !> \param blacs_env ...
173 ! **************************************************************************************************
174  SUBROUTINE cp_dbcsr_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa, &
175  para_env, blacs_env)
176  TYPE(dbcsr_type) :: matrix
177  INTEGER, INTENT(IN) :: neig
178  TYPE(dbcsr_type) :: matrixb, matrixout
179  CHARACTER(LEN=*), INTENT(IN) :: op
180  CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: pos, transa
181  TYPE(mp_para_env_type), POINTER :: para_env
182  TYPE(cp_blacs_env_type), POINTER :: blacs_env
183 
184  CHARACTER(len=*), PARAMETER :: routinen = 'cp_dbcsr_cholesky_restore'
185 
186  CHARACTER :: chol_pos, chol_transa
187  INTEGER :: handle, nfullcols_total, nfullrows_total
188  TYPE(cp_fm_struct_type), POINTER :: fm_struct
189  TYPE(cp_fm_type) :: fm_matrix, fm_matrixb, fm_matrixout
190 
191  CALL timeset(routinen, handle)
192 
193  NULLIFY (fm_struct)
194 
195  CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
196  CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
197  ncol_global=nfullcols_total, para_env=para_env)
198  CALL cp_fm_create(fm_matrix, fm_struct, name="fm_matrix")
199  CALL cp_fm_struct_release(fm_struct)
200 
201  CALL dbcsr_get_info(matrixb, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
202  CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
203  ncol_global=nfullcols_total, para_env=para_env)
204  CALL cp_fm_create(fm_matrixb, fm_struct, name="fm_matrixb")
205  CALL cp_fm_struct_release(fm_struct)
206 
207  CALL dbcsr_get_info(matrixout, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
208  CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
209  ncol_global=nfullcols_total, para_env=para_env)
210  CALL cp_fm_create(fm_matrixout, fm_struct, name="fm_matrixout")
211  CALL cp_fm_struct_release(fm_struct)
212 
213  CALL copy_dbcsr_to_fm(matrix, fm_matrix)
214  CALL copy_dbcsr_to_fm(matrixb, fm_matrixb)
215  !CALL copy_dbcsr_to_fm(matrixout, fm_matrixout)
216 
217  IF (op /= "SOLVE" .AND. op /= "MULTIPLY") &
218  cpabort("wrong argument op")
219 
220  IF (PRESENT(pos)) THEN
221  SELECT CASE (pos)
222  CASE ("LEFT")
223  chol_pos = 'L'
224  CASE ("RIGHT")
225  chol_pos = 'R'
226  CASE DEFAULT
227  cpabort("wrong argument pos")
228  END SELECT
229  ELSE
230  chol_pos = 'L'
231  END IF
232 
233  chol_transa = 'N'
234  IF (PRESENT(transa)) chol_transa = transa
235 
236  IF ((fm_matrix%use_sp .NEQV. fm_matrixb%use_sp) .OR. (fm_matrix%use_sp .NEQV. fm_matrixout%use_sp)) &
237  cpabort("not the same precision")
238 
239  CALL cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, chol_pos, chol_transa)
240 
241  CALL copy_fm_to_dbcsr(fm_matrixout, matrixout)
242 
243  CALL cp_fm_release(fm_matrix)
244  CALL cp_fm_release(fm_matrixb)
245  CALL cp_fm_release(fm_matrixout)
246 
247  CALL timestop(handle)
248 
249  END SUBROUTINE cp_dbcsr_cholesky_restore
250 
251 END MODULE cp_dbcsr_cholesky
252 
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
cp_dbcsr_cholesky
Interface to (sca)lapack for the Cholesky based procedures.
Definition: cp_dbcsr_cholesky.F:17
cp_dbcsr_cholesky::cp_dbcsr_cholesky_decompose
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition: cp_dbcsr_cholesky.F:61
cp_dbcsr_cholesky::cp_dbcsr_cholesky_restore
subroutine, public cp_dbcsr_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa, para_env, blacs_env)
...
Definition: cp_dbcsr_cholesky.F:176
cp_dbcsr_cholesky::cp_dbcsr_cholesky_invert
subroutine, public cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, upper_to_full)
used to replace the cholesky decomposition by the inverse
Definition: cp_dbcsr_cholesky.F:114
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition: cp_dbcsr_operations.F:118
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_potrf
subroutine, public cp_fm_potrf(fm_matrix, n)
Cholesky decomposition.
Definition: cp_fm_basic_linalg.F:3085
cp_fm_basic_linalg::cp_fm_upper_to_full
subroutine, public cp_fm_upper_to_full(matrix, work)
given an upper triangular matrix computes the corresponding full matrix
Definition: cp_fm_basic_linalg.F:1759
cp_fm_basic_linalg::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa)
...
Definition: cp_fm_basic_linalg.F:3162
cp_fm_basic_linalg::cp_fm_potri
subroutine, public cp_fm_potri(fm_matrix, n)
Invert trianguar matrix.
Definition: cp_fm_basic_linalg.F:3122
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23