Line data Source code
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 arnoldi iteration using dbcsr
10 : !> \par History
11 : !> 2014.09 created [Florian Schiffmann]
12 : !> \author Florian Schiffmann
13 : ! **************************************************************************************************
14 :
15 : MODULE arnoldi_api
16 : USE arnoldi_data_methods, ONLY: arnoldi_is_converged,&
17 : deallocate_arnoldi_data,&
18 : get_nrestart,&
19 : get_selected_ritz_val,&
20 : get_selected_ritz_vector,&
21 : select_evals,&
22 : set_arnoldi_initial_vector,&
23 : setup_arnoldi_data
24 : USE arnoldi_methods, ONLY: arnoldi_init,&
25 : arnoldi_iram,&
26 : build_subspace,&
27 : compute_evals,&
28 : gev_arnoldi_init,&
29 : gev_build_subspace,&
30 : gev_update_data
31 : USE arnoldi_types, ONLY: arnoldi_control_type,&
32 : arnoldi_data_type,&
33 : get_control,&
34 : m_x_v_vectors_type
35 : USE dbcsr_api, ONLY: &
36 : dbcsr_add, dbcsr_copy, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, &
37 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
38 : dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type
39 : USE dbcsr_vector, ONLY: create_col_vec_from_matrix,&
40 : create_replicated_col_vec_from_matrix,&
41 : create_replicated_row_vec_from_matrix,&
42 : dbcsr_matrix_colvec_multiply
43 : USE kinds, ONLY: dp
44 : USE message_passing, ONLY: mp_comm_type
45 : #include "../base/base_uses.f90"
46 :
47 : IMPLICIT NONE
48 :
49 : PRIVATE
50 :
51 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'arnoldi_api'
52 :
53 : PUBLIC :: arnoldi_ev, arnoldi_extremal, arnoldi_conjugate_gradient
54 : PUBLIC :: arnoldi_data_type, setup_arnoldi_data, deallocate_arnoldi_data
55 : PUBLIC :: set_arnoldi_initial_vector
56 : PUBLIC :: get_selected_ritz_val, get_selected_ritz_vector
57 :
58 : CONTAINS
59 :
60 : ! **************************************************************************************************
61 : !> \brief Driver routine for different arnoldi eigenvalue methods
62 : !> the selection which one is to be taken is made beforehand in the
63 : !> setup call passing the generalized_ev flag or not
64 : !> \param matrix ...
65 : !> \param arnoldi_data ...
66 : ! **************************************************************************************************
67 :
68 130277 : SUBROUTINE arnoldi_ev(matrix, arnoldi_data)
69 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
70 : TYPE(arnoldi_data_type) :: arnoldi_data
71 :
72 : TYPE(arnoldi_control_type), POINTER :: control
73 :
74 130277 : control => get_control(arnoldi_data)
75 :
76 130277 : IF (control%generalized_ev) THEN
77 3692 : CALL arnoldi_generalized_ev(matrix, arnoldi_data)
78 : ELSE
79 126585 : CALL arnoldi_normal_ev(matrix, arnoldi_data)
80 : END IF
81 :
82 130277 : END SUBROUTINE arnoldi_ev
83 :
84 : ! **************************************************************************************************
85 : !> \brief The main routine for arnoldi method to compute ritz values
86 : !> vectors of a matrix. Can take multiple matrices to solve
87 : !> ( M(N)*...*M(2)*M(1) )*v=v*e. A, B, ... have to be merged in a array of pointers
88 : !> arnoldi data has to be create with the setup routine and
89 : !> will contain on input all necessary information to start/restart
90 : !> the calculation. On output it contains all data
91 : !> \param matrix a pointer array to dbcsr_matrices. Multiplication order is as
92 : !> described above
93 : !> \param arnoldi_data On input data_type contains all information to start/restart
94 : !> an arnoldi iteration
95 : !> On output all data areas are filled to allow arbitrary post
96 : !> processing of the created subspace
97 : !> arnoldi_data has to be created with setup_arnoldi_data
98 : ! **************************************************************************************************
99 126585 : SUBROUTINE arnoldi_normal_ev(matrix, arnoldi_data)
100 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
101 : TYPE(arnoldi_data_type) :: arnoldi_data
102 :
103 : CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_normal_ev'
104 :
105 : INTEGER :: handle, i_loop, ncol_local, nrow_local
106 : TYPE(arnoldi_control_type), POINTER :: control
107 : TYPE(dbcsr_type), POINTER :: restart_vec
108 : TYPE(m_x_v_vectors_type) :: vectors
109 :
110 126585 : NULLIFY (restart_vec)
111 126585 : CALL timeset(routineN, handle)
112 :
113 : !prepare the vector like matrives needed in the matrix vector products, they will be reused throughout the iterations
114 126585 : CALL create_col_vec_from_matrix(vectors%input_vec, matrix(1)%matrix, 1)
115 126585 : CALL dbcsr_copy(vectors%result_vec, vectors%input_vec)
116 126585 : CALL create_replicated_col_vec_from_matrix(vectors%rep_col_vec, matrix(1)%matrix, 1)
117 126585 : CALL create_replicated_row_vec_from_matrix(vectors%rep_row_vec, matrix(1)%matrix, 1)
118 :
119 : ! Tells whether we have local data available on the processor (usually all in pcol 0 but even there can be some without data)
120 126585 : control => get_control(arnoldi_data)
121 126585 : CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
122 126585 : control%local_comp = ncol_local > 0 .AND. nrow_local > 0
123 :
124 127579 : DO i_loop = 0, get_nrestart(arnoldi_data)
125 :
126 126895 : IF (.NOT. control%iram .OR. i_loop == 0) THEN
127 : ! perform the standard arnoldi, if restarts are requested use the first (only makes sense if 1ev is requested)
128 126657 : IF (ASSOCIATED(restart_vec)) CALL set_arnoldi_initial_vector(arnoldi_data, restart_vec)
129 126657 : CALL arnoldi_init(matrix, vectors, arnoldi_data)
130 : ELSE
131 : ! perform an implicit restart
132 238 : CALL arnoldi_iram(arnoldi_data)
133 : END IF
134 :
135 : ! Generate the subspace
136 126895 : CALL build_subspace(matrix, vectors, arnoldi_data)
137 :
138 : ! If we reached the maximum number of steps or the subspace converged we still need to get the eigenvalues
139 126895 : CALL compute_evals(arnoldi_data)
140 :
141 : ! Select the evals according to user selection and keep them in arnoldi_data
142 126895 : CALL select_evals(arnoldi_data)
143 :
144 : ! Prepare for a restart with the best eigenvector not needed in case of iram but who cares
145 126895 : IF (.NOT. ASSOCIATED(restart_vec)) ALLOCATE (restart_vec)
146 126895 : CALL get_selected_ritz_vector(arnoldi_data, 1, matrix(1)%matrix, restart_vec)
147 :
148 : ! Check whether we can already go home
149 127579 : IF (control%converged) EXIT
150 : END DO
151 :
152 : ! Deallocated the work vectors
153 126585 : CALL dbcsr_release(vectors%input_vec)
154 126585 : CALL dbcsr_release(vectors%result_vec)
155 126585 : CALL dbcsr_release(vectors%rep_col_vec)
156 126585 : CALL dbcsr_release(vectors%rep_row_vec)
157 126585 : CALL dbcsr_release(restart_vec)
158 126585 : DEALLOCATE (restart_vec)
159 126585 : CALL timestop(handle)
160 :
161 126585 : END SUBROUTINE arnoldi_normal_ev
162 :
163 : ! **************************************************************************************************
164 : !> \brief The main routine for arnoldi method to compute the lowest ritz pair
165 : !> of a symmetric generalized eigenvalue problem.
166 : !> as input it takes a vector of matrices which for the GEV:
167 : !> M(1)*v=M(2)*v*lambda
168 : !> In other words, M(1) is the matrix and M(2) the metric
169 : !> This only works if the two matrices are symmetric in values
170 : !> (flag in dbcsr does not need to be set)
171 : !> \param matrix a pointer array to dbcsr_matrices. Multiplication order is as
172 : !> described above
173 : !> \param arnoldi_data On input data_type contains all information to start/restart
174 : !> an arnoldi iteration
175 : !> On output all data areas are filled to allow arbitrary post
176 : !> processing of the created subspace
177 : !> arnoldi_data has to be created with setup_arnoldi_data
178 : ! **************************************************************************************************
179 3692 : SUBROUTINE arnoldi_generalized_ev(matrix, arnoldi_data)
180 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
181 : TYPE(arnoldi_data_type) :: arnoldi_data
182 :
183 : CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_generalized_ev'
184 :
185 : INTEGER :: handle, i_loop, ncol_local, nrow_local
186 : TYPE(arnoldi_control_type), POINTER :: control
187 3692 : TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:) :: matrix_arnoldi
188 : TYPE(dbcsr_type), TARGET :: A_rho_B
189 : TYPE(m_x_v_vectors_type) :: vectors
190 :
191 3692 : CALL timeset(routineN, handle)
192 11076 : ALLOCATE (matrix_arnoldi(2))
193 : ! this matrix will contain +/- A-rho*B
194 3692 : matrix_arnoldi(1)%matrix => A_rho_B
195 3692 : matrix_arnoldi(2)%matrix => matrix(2)%matrix
196 :
197 : !prepare the vector like matrives needed in the matrix vector products, they will be reused throughout the iterations
198 3692 : CALL create_col_vec_from_matrix(vectors%input_vec, matrix(1)%matrix, 1)
199 3692 : CALL dbcsr_copy(vectors%result_vec, vectors%input_vec)
200 3692 : CALL create_replicated_col_vec_from_matrix(vectors%rep_col_vec, matrix(1)%matrix, 1)
201 3692 : CALL create_replicated_row_vec_from_matrix(vectors%rep_row_vec, matrix(1)%matrix, 1)
202 :
203 : ! Tells whether we have local data available on the processor (usually all in pcol 0 but even there can be some without data)
204 3692 : control => get_control(arnoldi_data)
205 3692 : CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
206 3692 : control%local_comp = ncol_local > 0 .AND. nrow_local > 0
207 :
208 4567 : DO i_loop = 0, get_nrestart(arnoldi_data)
209 4567 : IF (i_loop == 0) THEN
210 : ! perform the standard arnoldi initialization with a random vector
211 3692 : CALL gev_arnoldi_init(matrix, matrix_arnoldi, vectors, arnoldi_data)
212 : END IF
213 :
214 : ! Generate the subspace
215 4567 : CALL gev_build_subspace(matrix_arnoldi, vectors, arnoldi_data)
216 :
217 : ! If we reached the maximum number of steps or the subspace converged we still need to get the eigenvalues
218 4567 : CALL compute_evals(arnoldi_data)
219 :
220 : ! Select the evals according to user selection and keep them in arnoldi_data
221 4567 : CALL select_evals(arnoldi_data)
222 :
223 : ! update the matrices and compute the convergence
224 4567 : CALL gev_update_data(matrix, matrix_arnoldi, vectors, arnoldi_data)
225 :
226 : ! Check whether we can already go home
227 4567 : IF (control%converged) EXIT
228 : END DO
229 :
230 : ! Deallocated the work vectors
231 3692 : CALL dbcsr_release(vectors%input_vec)
232 3692 : CALL dbcsr_release(vectors%result_vec)
233 3692 : CALL dbcsr_release(vectors%rep_col_vec)
234 3692 : CALL dbcsr_release(vectors%rep_row_vec)
235 3692 : CALL dbcsr_release(A_rho_B)
236 3692 : DEALLOCATE (matrix_arnoldi)
237 :
238 3692 : CALL timestop(handle)
239 :
240 11076 : END SUBROUTINE arnoldi_generalized_ev
241 :
242 : ! **************************************************************************************************
243 : !> \brief simple wrapper to estimate extremal eigenvalues with arnoldi, using the old lanczos interface
244 : !> this hides some of the power of the arnoldi routines (e.g. only min or max eval, generalized eval, ritz vectors, etc.),
245 : !> and does not allow for providing an initial guess of the ritz vector.
246 : !> \param matrix_a input mat
247 : !> \param max_ev estimated max eval
248 : !> \param min_ev estimated min eval
249 : !> \param converged Usually arnoldi is more accurate than claimed.
250 : !> \param threshold target precision
251 : !> \param max_iter max allowed iterations (will be rounded up)
252 : ! **************************************************************************************************
253 122195 : SUBROUTINE arnoldi_extremal(matrix_a, max_ev, min_ev, converged, threshold, max_iter)
254 : TYPE(dbcsr_type), INTENT(INOUT), TARGET :: matrix_a
255 : REAL(KIND=dp), INTENT(OUT) :: max_ev, min_ev
256 : LOGICAL, INTENT(OUT) :: converged
257 : REAL(KIND=dp), INTENT(IN) :: threshold
258 : INTEGER, INTENT(IN) :: max_iter
259 :
260 : CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_extremal'
261 :
262 : INTEGER :: handle, max_iter_internal, nrestarts
263 : TYPE(arnoldi_data_type) :: my_arnoldi
264 122195 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: arnoldi_matrices
265 :
266 122195 : CALL timeset(routineN, handle)
267 :
268 : ! we go in chunks of max_iter_internal, and restart ater each of those.
269 : ! at low threshold smaller values of max_iter_internal make sense
270 122195 : IF (.TRUE.) max_iter_internal = 16
271 122195 : IF (threshold <= 1.0E-3_dp) max_iter_internal = 32
272 122195 : IF (threshold <= 1.0E-4_dp) max_iter_internal = 64
273 :
274 : ! the max number of iter will be (nrestarts+1)*max_iter_internal
275 122195 : nrestarts = max_iter/max_iter_internal
276 :
277 244390 : ALLOCATE (arnoldi_matrices(1))
278 122195 : arnoldi_matrices(1)%matrix => matrix_a
279 : CALL setup_arnoldi_data(my_arnoldi, arnoldi_matrices, max_iter=max_iter_internal, &
280 : threshold=threshold, selection_crit=1, nval_request=2, nrestarts=nrestarts, &
281 122195 : generalized_ev=.FALSE., iram=.TRUE.)
282 122195 : CALL arnoldi_ev(arnoldi_matrices, my_arnoldi)
283 122195 : converged = arnoldi_is_converged(my_arnoldi)
284 122195 : max_eV = REAL(get_selected_ritz_val(my_arnoldi, 2), dp)
285 122195 : min_eV = REAL(get_selected_ritz_val(my_arnoldi, 1), dp)
286 122195 : CALL deallocate_arnoldi_data(my_arnoldi)
287 122195 : DEALLOCATE (arnoldi_matrices)
288 :
289 122195 : CALL timestop(handle)
290 :
291 122195 : END SUBROUTINE arnoldi_extremal
292 :
293 : ! **************************************************************************************************
294 : !> \brief Wrapper for conjugated gradient algorithm for Ax=b
295 : !> \param matrix_a input mat
296 : !> \param vec_x input right hand side vector; output solution vector, fully replicated!
297 : !> \param matrix_p input preconditioner (optional)
298 : !> \param converged ...
299 : !> \param threshold target precision
300 : !> \param max_iter max allowed iterations
301 : ! **************************************************************************************************
302 198 : SUBROUTINE arnoldi_conjugate_gradient(matrix_a, vec_x, matrix_p, converged, threshold, max_iter)
303 : TYPE(dbcsr_type), INTENT(IN), TARGET :: matrix_a
304 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: vec_x
305 : TYPE(dbcsr_type), INTENT(IN), OPTIONAL, TARGET :: matrix_p
306 : LOGICAL, INTENT(OUT) :: converged
307 : REAL(KIND=dp), INTENT(IN) :: threshold
308 : INTEGER, INTENT(IN) :: max_iter
309 :
310 : CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_conjugate_gradient'
311 :
312 : INTEGER :: handle, i, j, nb, nloc, no
313 198 : INTEGER, DIMENSION(:), POINTER :: rb_offset, rb_size
314 198 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xvec
315 : TYPE(arnoldi_control_type), POINTER :: control
316 : TYPE(arnoldi_data_type) :: my_arnoldi
317 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
318 198 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: arnoldi_matrices
319 : TYPE(dbcsr_type), TARGET :: x
320 : TYPE(m_x_v_vectors_type) :: vectors
321 :
322 198 : CALL timeset(routineN, handle)
323 :
324 : !prepare the vector like matrices needed in the matrix vector products,
325 : !they will be reused throughout the iterations
326 198 : CALL create_col_vec_from_matrix(vectors%input_vec, matrix_a, 1)
327 198 : CALL dbcsr_copy(vectors%result_vec, vectors%input_vec)
328 198 : CALL create_replicated_col_vec_from_matrix(vectors%rep_col_vec, matrix_a, 1)
329 198 : CALL create_replicated_row_vec_from_matrix(vectors%rep_row_vec, matrix_a, 1)
330 : !
331 198 : CALL dbcsr_copy(x, vectors%input_vec)
332 : !
333 198 : CALL dbcsr_get_info(x, nfullrows_local=nloc, row_blk_size=rb_size, row_blk_offset=rb_offset)
334 : !
335 198 : CALL dbcsr_iterator_start(dbcsr_iter, x)
336 607 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
337 409 : CALL dbcsr_iterator_next_block(dbcsr_iter, i, j, xvec)
338 409 : nb = rb_size(i)
339 409 : no = rb_offset(i)
340 4984 : xvec(1:nb, 1) = vec_x(no:no + nb - 1)
341 : END DO
342 198 : CALL dbcsr_iterator_stop(dbcsr_iter)
343 :
344 : ! Arnoldi interface (used just for the iterator information here
345 792 : ALLOCATE (arnoldi_matrices(3))
346 198 : arnoldi_matrices(1)%matrix => matrix_a
347 198 : IF (PRESENT(matrix_p)) THEN
348 198 : arnoldi_matrices(2)%matrix => matrix_p
349 : ELSE
350 0 : NULLIFY (arnoldi_matrices(2)%matrix)
351 : END IF
352 198 : arnoldi_matrices(3)%matrix => x
353 : CALL setup_arnoldi_data(my_arnoldi, arnoldi_matrices, max_iter=max_iter, &
354 : threshold=threshold, selection_crit=1, nval_request=1, nrestarts=0, &
355 198 : generalized_ev=.FALSE., iram=.FALSE.)
356 :
357 198 : CALL conjugate_gradient(my_arnoldi, arnoldi_matrices, vectors)
358 :
359 198 : converged = arnoldi_is_converged(my_arnoldi)
360 :
361 8952 : vec_x = 0.0_dp
362 198 : CALL dbcsr_iterator_start(dbcsr_iter, x)
363 607 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
364 409 : CALL dbcsr_iterator_next_block(dbcsr_iter, i, j, xvec)
365 409 : nb = rb_size(i)
366 409 : no = rb_offset(i)
367 4984 : vec_x(no:no + nb - 1) = xvec(1:nb, 1)
368 : END DO
369 198 : CALL dbcsr_iterator_stop(dbcsr_iter)
370 198 : control => get_control(my_arnoldi)
371 17706 : CALL control%mp_group%sum(vec_x)
372 : ! Deallocated the work vectors
373 198 : CALL dbcsr_release(x)
374 198 : CALL dbcsr_release(vectors%input_vec)
375 198 : CALL dbcsr_release(vectors%result_vec)
376 198 : CALL dbcsr_release(vectors%rep_col_vec)
377 198 : CALL dbcsr_release(vectors%rep_row_vec)
378 :
379 198 : CALL deallocate_arnoldi_data(my_arnoldi)
380 198 : DEALLOCATE (arnoldi_matrices)
381 :
382 198 : CALL timestop(handle)
383 :
384 396 : END SUBROUTINE arnoldi_conjugate_gradient
385 :
386 : ! **************************************************************************************************
387 : !> \brief ...
388 : !> \param arnoldi_data ...
389 : !> \param arnoldi_matrices ...
390 : !> \param vectors ...
391 : ! **************************************************************************************************
392 198 : SUBROUTINE conjugate_gradient(arnoldi_data, arnoldi_matrices, vectors)
393 : TYPE(arnoldi_data_type) :: arnoldi_data
394 : TYPE(dbcsr_p_type), DIMENSION(:) :: arnoldi_matrices
395 : TYPE(m_x_v_vectors_type) :: vectors
396 :
397 : INTEGER :: iter
398 : REAL(KIND=dp) :: alpha, beta, pap, rsnew, rsold
399 : TYPE(arnoldi_control_type), POINTER :: control
400 : TYPE(dbcsr_type) :: apvec, pvec, rvec, zvec
401 : TYPE(dbcsr_type), POINTER :: amat, pmat, xvec
402 : TYPE(mp_comm_type) :: mpgrp, pcgrp
403 :
404 396 : control => get_control(arnoldi_data)
405 198 : control%converged = .FALSE.
406 198 : pcgrp = control%pcol_group
407 198 : mpgrp = control%mp_group
408 :
409 198 : NULLIFY (amat, pmat, xvec)
410 198 : amat => arnoldi_matrices(1)%matrix
411 198 : pmat => arnoldi_matrices(2)%matrix
412 198 : xvec => arnoldi_matrices(3)%matrix
413 :
414 198 : IF (ASSOCIATED(pmat)) THEN
415 : ! Preconditioned conjugate gradients
416 198 : CALL dbcsr_copy(vectors%input_vec, xvec)
417 : CALL dbcsr_matrix_colvec_multiply(pmat, vectors%input_vec, vectors%result_vec, 1.0_dp, &
418 198 : 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec)
419 198 : CALL dbcsr_copy(pvec, vectors%result_vec)
420 198 : CALL dbcsr_copy(vectors%input_vec, pvec)
421 : CALL dbcsr_matrix_colvec_multiply(amat, vectors%input_vec, vectors%result_vec, 1.0_dp, &
422 198 : 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec)
423 198 : CALL dbcsr_copy(apvec, vectors%result_vec)
424 198 : CALL dbcsr_copy(rvec, xvec)
425 198 : CALL dbcsr_add(rvec, apvec, 1.0_dp, -1.0_dp)
426 198 : CALL dbcsr_copy(xvec, pvec)
427 198 : CALL dbcsr_copy(vectors%input_vec, rvec)
428 : CALL dbcsr_matrix_colvec_multiply(pmat, vectors%input_vec, vectors%result_vec, 1.0_dp, &
429 198 : 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec)
430 198 : CALL dbcsr_copy(zvec, vectors%result_vec)
431 198 : CALL dbcsr_copy(pvec, zvec)
432 198 : rsold = vec_dot_vec(rvec, zvec, mpgrp)
433 5350 : DO iter = 1, control%max_iter
434 5338 : CALL dbcsr_copy(vectors%input_vec, pvec)
435 : CALL dbcsr_matrix_colvec_multiply(amat, vectors%input_vec, vectors%result_vec, 1.0_dp, &
436 5338 : 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec)
437 5338 : CALL dbcsr_copy(apvec, vectors%result_vec)
438 :
439 5338 : pap = vec_dot_vec(pvec, apvec, mpgrp)
440 5338 : IF (ABS(pap) < 1.e-24_dp) THEN
441 18 : alpha = 0.0_dp
442 : ELSE
443 5320 : alpha = rsold/pap
444 : END IF
445 :
446 5338 : CALL dbcsr_add(xvec, pvec, 1.0_dp, alpha)
447 5338 : CALL dbcsr_add(rvec, apvec, 1.0_dp, -alpha)
448 5338 : rsnew = vec_dot_vec(rvec, rvec, mpgrp)
449 5338 : IF (SQRT(rsnew) < control%threshold) EXIT
450 5152 : CPASSERT(alpha /= 0.0_dp)
451 :
452 5152 : CALL dbcsr_copy(vectors%input_vec, rvec)
453 : CALL dbcsr_matrix_colvec_multiply(pmat, vectors%input_vec, vectors%result_vec, 1.0_dp, &
454 5152 : 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec)
455 5152 : CALL dbcsr_copy(zvec, vectors%result_vec)
456 5152 : rsnew = vec_dot_vec(rvec, zvec, mpgrp)
457 5152 : beta = rsnew/rsold
458 5152 : CALL dbcsr_add(pvec, zvec, beta, 1.0_dp)
459 5350 : rsold = rsnew
460 : END DO
461 198 : IF (SQRT(rsnew) < control%threshold) control%converged = .TRUE.
462 198 : CALL dbcsr_release(zvec)
463 198 : CALL dbcsr_release(pvec)
464 198 : CALL dbcsr_release(rvec)
465 198 : CALL dbcsr_release(apvec)
466 :
467 : ELSE
468 0 : CALL dbcsr_copy(pvec, xvec)
469 0 : CALL dbcsr_copy(rvec, xvec)
470 0 : CALL dbcsr_set(xvec, 0.0_dp)
471 : ! Conjugate gradients
472 0 : rsold = vec_dot_vec(rvec, rvec, mpgrp)
473 0 : DO iter = 1, control%max_iter
474 0 : CALL dbcsr_copy(vectors%input_vec, pvec)
475 : CALL dbcsr_matrix_colvec_multiply(amat, vectors%input_vec, vectors%result_vec, 1.0_dp, &
476 0 : 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec)
477 0 : CALL dbcsr_copy(apvec, vectors%result_vec)
478 0 : pap = vec_dot_vec(pvec, apvec, mpgrp)
479 0 : IF (ABS(pap) < 1.e-24_dp) THEN
480 0 : alpha = 0.0_dp
481 : ELSE
482 0 : alpha = rsold/pap
483 : END IF
484 0 : CALL dbcsr_add(xvec, pvec, 1.0_dp, alpha)
485 0 : CALL dbcsr_add(rvec, apvec, 1.0_dp, -alpha)
486 0 : rsnew = vec_dot_vec(rvec, rvec, mpgrp)
487 0 : IF (SQRT(rsnew) < control%threshold) EXIT
488 0 : CPASSERT(alpha /= 0.0_dp)
489 0 : beta = rsnew/rsold
490 0 : CALL dbcsr_add(pvec, rvec, beta, 1.0_dp)
491 0 : rsold = rsnew
492 : END DO
493 0 : IF (SQRT(rsnew) < control%threshold) control%converged = .TRUE.
494 0 : CALL dbcsr_release(pvec)
495 0 : CALL dbcsr_release(rvec)
496 0 : CALL dbcsr_release(apvec)
497 : END IF
498 :
499 198 : END SUBROUTINE conjugate_gradient
500 :
501 : ! **************************************************************************************************
502 : !> \brief ...
503 : !> \param avec ...
504 : !> \param bvec ...
505 : !> \param mpgrp ...
506 : !> \return ...
507 : ! **************************************************************************************************
508 16026 : FUNCTION vec_dot_vec(avec, bvec, mpgrp) RESULT(adotb)
509 : TYPE(dbcsr_type) :: avec, bvec
510 : TYPE(mp_comm_type), INTENT(IN) :: mpgrp
511 : REAL(KIND=dp) :: adotb
512 :
513 : INTEGER :: i, j
514 : LOGICAL :: found
515 16026 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: av, bv
516 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
517 :
518 16026 : adotb = 0.0_dp
519 16026 : CALL dbcsr_iterator_start(dbcsr_iter, avec)
520 66930 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
521 50904 : CALL dbcsr_iterator_next_block(dbcsr_iter, i, j, av)
522 50904 : CALL dbcsr_get_block_p(bvec, i, j, bv, found)
523 66930 : IF (found .AND. SIZE(bv) > 0) THEN
524 955968 : adotb = adotb + DOT_PRODUCT(av(:, 1), bv(:, 1))
525 : END IF
526 : END DO
527 16026 : CALL dbcsr_iterator_stop(dbcsr_iter)
528 16026 : CALL mpgrp%sum(adotb)
529 :
530 16026 : END FUNCTION vec_dot_vec
531 :
532 : END MODULE arnoldi_api
|
