(git:ed6f26b)
Loading...
Searching...
No Matches
arnoldi_api.F
Go to the documentation of this file.
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 arnoldi iteration using dbcsr
10!> \par History
11!> 2014.09 created [Florian Schiffmann]
12!> 2023.12 Removed support for single-precision [Ole Schuett]
13!> 2024.12 Removed support for complex input matrices [Ole Schuett]
14!> \author Florian Schiffmann
15! **************************************************************************************************
16MODULE arnoldi_api
17 USE arnoldi_data_methods, ONLY: arnoldi_is_converged,&
18 deallocate_arnoldi_env,&
19 get_nrestart,&
20 get_selected_ritz_val,&
21 get_selected_ritz_vector,&
22 select_evals,&
23 set_arnoldi_initial_vector,&
24 setup_arnoldi_env
25 USE arnoldi_methods, ONLY: arnoldi_init,&
26 arnoldi_iram,&
27 build_subspace,&
28 compute_evals,&
29 gev_arnoldi_init,&
30 gev_build_subspace,&
31 gev_update_data
32 USE arnoldi_types, ONLY: arnoldi_control_type,&
33 arnoldi_env_type,&
34 get_control,&
35 m_x_v_vectors_type
36 USE arnoldi_vector, ONLY: create_col_vec_from_matrix,&
37 create_replicated_col_vec_from_matrix,&
38 create_replicated_row_vec_from_matrix,&
39 dbcsr_matrix_colvec_multiply
40 USE cp_dbcsr_api, ONLY: &
41 dbcsr_add, dbcsr_copy, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, &
42 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
43 dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type
44 USE kinds, ONLY: dp
45 USE message_passing, ONLY: mp_comm_type
46#include "../base/base_uses.f90"
47
48 IMPLICIT NONE
49
50 PRIVATE
51
52 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'arnoldi_api'
53
54 PUBLIC :: arnoldi_ev, arnoldi_extremal, arnoldi_conjugate_gradient
55 PUBLIC :: arnoldi_env_type, setup_arnoldi_env, deallocate_arnoldi_env
56 PUBLIC :: set_arnoldi_initial_vector
57 PUBLIC :: get_selected_ritz_val, get_selected_ritz_vector
58
59CONTAINS
60
61! **************************************************************************************************
62!> \brief Driver routine for different arnoldi eigenvalue methods
63!> the selection which one is to be taken is made beforehand in the
64!> setup call passing the generalized_ev flag or not
65!> \param matrix ...
66!> \param arnoldi_env ...
67! **************************************************************************************************
68 SUBROUTINE arnoldi_ev(matrix, arnoldi_env)
69 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
70 TYPE(arnoldi_env_type) :: arnoldi_env
71
72 TYPE(arnoldi_control_type), POINTER :: control
73
74 control => get_control(arnoldi_env)
75
76 IF (control%generalized_ev) THEN
77 CALL arnoldi_generalized_ev(matrix, arnoldi_env)
78 ELSE
79 CALL arnoldi_normal_ev(matrix, arnoldi_env)
80 END IF
81
82 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_env 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_env has to be created with setup_arnoldi_env
98! **************************************************************************************************
99 SUBROUTINE arnoldi_normal_ev(matrix, arnoldi_env)
100 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
101 TYPE(arnoldi_env_type) :: arnoldi_env
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 NULLIFY (restart_vec)
111 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 CALL create_col_vec_from_matrix(vectors%input_vec, matrix(1)%matrix, 1)
115 CALL dbcsr_copy(vectors%result_vec, vectors%input_vec)
116 CALL create_replicated_col_vec_from_matrix(vectors%rep_col_vec, matrix(1)%matrix, 1)
117 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 control => get_control(arnoldi_env)
121 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
122 control%local_comp = ncol_local > 0 .AND. nrow_local > 0
123
124 DO i_loop = 0, get_nrestart(arnoldi_env)
125
126 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 IF (ASSOCIATED(restart_vec)) CALL set_arnoldi_initial_vector(arnoldi_env, restart_vec)
129 CALL arnoldi_init(matrix, vectors, arnoldi_env)
130 ELSE
131 ! perform an implicit restart
132 CALL arnoldi_iram(arnoldi_env)
133 END IF
134
135 ! Generate the subspace
136 CALL build_subspace(matrix, vectors, arnoldi_env)
137
138 ! If we reached the maximum number of steps or the subspace converged we still need to get the eigenvalues
139 CALL compute_evals(arnoldi_env)
140
141 ! Select the evals according to user selection and keep them in arnoldi_env
142 CALL select_evals(arnoldi_env)
143
144 ! Prepare for a restart with the best eigenvector not needed in case of iram but who cares
145 IF (.NOT. ASSOCIATED(restart_vec)) ALLOCATE (restart_vec)
146 CALL get_selected_ritz_vector(arnoldi_env, 1, matrix(1)%matrix, restart_vec)
147
148 ! Check whether we can already go home
149 IF (control%converged) EXIT
150 END DO
151
152 ! Deallocated the work vectors
153 CALL dbcsr_release(vectors%input_vec)
154 CALL dbcsr_release(vectors%result_vec)
155 CALL dbcsr_release(vectors%rep_col_vec)
156 CALL dbcsr_release(vectors%rep_row_vec)
157 CALL dbcsr_release(restart_vec)
158 DEALLOCATE (restart_vec)
159 CALL timestop(handle)
160
161 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_env 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_env has to be created with setup_arnoldi_env
178! **************************************************************************************************
179 SUBROUTINE arnoldi_generalized_ev(matrix, arnoldi_env)
180 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
181 TYPE(arnoldi_env_type) :: arnoldi_env
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 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 CALL timeset(routinen, handle)
192 ALLOCATE (matrix_arnoldi(2))
193 ! this matrix will contain +/- A-rho*B
194 matrix_arnoldi(1)%matrix => a_rho_b
195 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 CALL create_col_vec_from_matrix(vectors%input_vec, matrix(1)%matrix, 1)
199 CALL dbcsr_copy(vectors%result_vec, vectors%input_vec)
200 CALL create_replicated_col_vec_from_matrix(vectors%rep_col_vec, matrix(1)%matrix, 1)
201 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 control => get_control(arnoldi_env)
205 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
206 control%local_comp = ncol_local > 0 .AND. nrow_local > 0
207
208 DO i_loop = 0, get_nrestart(arnoldi_env)
209 IF (i_loop == 0) THEN
210 ! perform the standard arnoldi initialization with a random vector
211 CALL gev_arnoldi_init(matrix, matrix_arnoldi, vectors, arnoldi_env)
212 END IF
213
214 ! Generate the subspace
215 CALL gev_build_subspace(matrix_arnoldi, vectors, arnoldi_env)
216
217 ! If we reached the maximum number of steps or the subspace converged we still need to get the eigenvalues
218 CALL compute_evals(arnoldi_env)
219
220 ! Select the evals according to user selection and keep them in arnoldi_env
221 CALL select_evals(arnoldi_env)
222
223 ! update the matrices and compute the convergence
224 CALL gev_update_data(matrix, matrix_arnoldi, vectors, arnoldi_env)
225
226 ! Check whether we can already go home
227 IF (control%converged) EXIT
228 END DO
229
230 ! Deallocated the work vectors
231 CALL dbcsr_release(vectors%input_vec)
232 CALL dbcsr_release(vectors%result_vec)
233 CALL dbcsr_release(vectors%rep_col_vec)
234 CALL dbcsr_release(vectors%rep_row_vec)
235 CALL dbcsr_release(a_rho_b)
236 DEALLOCATE (matrix_arnoldi)
237
238 CALL timestop(handle)
239
240 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 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_env_type) :: arnoldi_env
264 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: arnoldi_matrices
265
266 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 IF (.true.) max_iter_internal = 16
271 IF (threshold <= 1.0e-3_dp) max_iter_internal = 32
272 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 nrestarts = max_iter/max_iter_internal
276
277 ALLOCATE (arnoldi_matrices(1))
278 arnoldi_matrices(1)%matrix => matrix_a
279 CALL setup_arnoldi_env(arnoldi_env, arnoldi_matrices, max_iter=max_iter_internal, &
280 threshold=threshold, selection_crit=1, nval_request=2, nrestarts=nrestarts, &
281 generalized_ev=.false., iram=.true.)
282 CALL arnoldi_ev(arnoldi_matrices, arnoldi_env)
283 converged = arnoldi_is_converged(arnoldi_env)
284 max_ev = real(get_selected_ritz_val(arnoldi_env, 2), dp)
285 min_ev = real(get_selected_ritz_val(arnoldi_env, 1), dp)
286 CALL deallocate_arnoldi_env(arnoldi_env)
287 DEALLOCATE (arnoldi_matrices)
288
289 CALL timestop(handle)
290
291 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 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 INTEGER, DIMENSION(:), POINTER :: rb_offset, rb_size
314 REAL(kind=dp), DIMENSION(:, :), POINTER :: xvec
315 TYPE(arnoldi_control_type), POINTER :: control
316 TYPE(arnoldi_env_type) :: arnoldi_env
317 TYPE(dbcsr_iterator_type) :: dbcsr_iter
318 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: arnoldi_matrices
319 TYPE(dbcsr_type), TARGET :: x
320 TYPE(m_x_v_vectors_type) :: vectors
321
322 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 CALL create_col_vec_from_matrix(vectors%input_vec, matrix_a, 1)
327 CALL dbcsr_copy(vectors%result_vec, vectors%input_vec)
328 CALL create_replicated_col_vec_from_matrix(vectors%rep_col_vec, matrix_a, 1)
329 CALL create_replicated_row_vec_from_matrix(vectors%rep_row_vec, matrix_a, 1)
330 !
331 CALL dbcsr_copy(x, vectors%input_vec)
332 !
333 CALL dbcsr_get_info(x, nfullrows_local=nloc, row_blk_size=rb_size, row_blk_offset=rb_offset)
334 !
335 CALL dbcsr_iterator_start(dbcsr_iter, x)
336 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
337 CALL dbcsr_iterator_next_block(dbcsr_iter, i, j, xvec)
338 nb = rb_size(i)
339 no = rb_offset(i)
340 xvec(1:nb, 1) = vec_x(no:no + nb - 1)
341 END DO
342 CALL dbcsr_iterator_stop(dbcsr_iter)
343
344 ! Arnoldi interface (used just for the iterator information here
345 ALLOCATE (arnoldi_matrices(3))
346 arnoldi_matrices(1)%matrix => matrix_a
347 IF (PRESENT(matrix_p)) THEN
348 arnoldi_matrices(2)%matrix => matrix_p
349 ELSE
350 NULLIFY (arnoldi_matrices(2)%matrix)
351 END IF
352 arnoldi_matrices(3)%matrix => x
353 CALL setup_arnoldi_env(arnoldi_env, arnoldi_matrices, max_iter=max_iter, &
354 threshold=threshold, selection_crit=1, nval_request=1, nrestarts=0, &
355 generalized_ev=.false., iram=.false.)
356
357 CALL conjugate_gradient(arnoldi_env, arnoldi_matrices, vectors)
358
359 converged = arnoldi_is_converged(arnoldi_env)
360
361 vec_x = 0.0_dp
362 CALL dbcsr_iterator_start(dbcsr_iter, x)
363 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
364 CALL dbcsr_iterator_next_block(dbcsr_iter, i, j, xvec)
365 nb = rb_size(i)
366 no = rb_offset(i)
367 vec_x(no:no + nb - 1) = xvec(1:nb, 1)
368 END DO
369 CALL dbcsr_iterator_stop(dbcsr_iter)
370 control => get_control(arnoldi_env)
371 CALL control%mp_group%sum(vec_x)
372 ! Deallocated the work vectors
373 CALL dbcsr_release(x)
374 CALL dbcsr_release(vectors%input_vec)
375 CALL dbcsr_release(vectors%result_vec)
376 CALL dbcsr_release(vectors%rep_col_vec)
377 CALL dbcsr_release(vectors%rep_row_vec)
378
379 CALL deallocate_arnoldi_env(arnoldi_env)
380 DEALLOCATE (arnoldi_matrices)
381
382 CALL timestop(handle)
383
384 END SUBROUTINE arnoldi_conjugate_gradient
385
386! **************************************************************************************************
387!> \brief ...
388!> \param arnoldi_env ...
389!> \param arnoldi_matrices ...
390!> \param vectors ...
391! **************************************************************************************************
392 SUBROUTINE conjugate_gradient(arnoldi_env, arnoldi_matrices, vectors)
393 TYPE(arnoldi_env_type) :: arnoldi_env
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 control => get_control(arnoldi_env)
405 control%converged = .false.
406 pcgrp = control%pcol_group
407 mpgrp = control%mp_group
408
409 NULLIFY (amat, pmat, xvec)
410 amat => arnoldi_matrices(1)%matrix
411 pmat => arnoldi_matrices(2)%matrix
412 xvec => arnoldi_matrices(3)%matrix
413
414 IF (ASSOCIATED(pmat)) THEN
415 ! Preconditioned conjugate gradients
416 CALL dbcsr_copy(vectors%input_vec, xvec)
417 CALL dbcsr_matrix_colvec_multiply(pmat, vectors%input_vec, vectors%result_vec, 1.0_dp, &
418 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec)
419 CALL dbcsr_copy(pvec, vectors%result_vec)
420 CALL dbcsr_copy(vectors%input_vec, pvec)
421 CALL dbcsr_matrix_colvec_multiply(amat, vectors%input_vec, vectors%result_vec, 1.0_dp, &
422 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec)
423 CALL dbcsr_copy(apvec, vectors%result_vec)
424 CALL dbcsr_copy(rvec, xvec)
425 CALL dbcsr_add(rvec, apvec, 1.0_dp, -1.0_dp)
426 CALL dbcsr_copy(xvec, pvec)
427 CALL dbcsr_copy(vectors%input_vec, rvec)
428 CALL dbcsr_matrix_colvec_multiply(pmat, vectors%input_vec, vectors%result_vec, 1.0_dp, &
429 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec)
430 CALL dbcsr_copy(zvec, vectors%result_vec)
431 CALL dbcsr_copy(pvec, zvec)
432 rsold = vec_dot_vec(rvec, zvec, mpgrp)
433 DO iter = 1, control%max_iter
434 CALL dbcsr_copy(vectors%input_vec, pvec)
435 CALL dbcsr_matrix_colvec_multiply(amat, vectors%input_vec, vectors%result_vec, 1.0_dp, &
436 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec)
437 CALL dbcsr_copy(apvec, vectors%result_vec)
438
439 pap = vec_dot_vec(pvec, apvec, mpgrp)
440 IF (abs(pap) < 1.e-24_dp) THEN
441 alpha = 0.0_dp
442 ELSE
443 alpha = rsold/pap
444 END IF
445
446 CALL dbcsr_add(xvec, pvec, 1.0_dp, alpha)
447 CALL dbcsr_add(rvec, apvec, 1.0_dp, -alpha)
448 rsnew = vec_dot_vec(rvec, rvec, mpgrp)
449 IF (sqrt(rsnew) < control%threshold) EXIT
450 cpassert(alpha /= 0.0_dp)
451
452 CALL dbcsr_copy(vectors%input_vec, rvec)
453 CALL dbcsr_matrix_colvec_multiply(pmat, vectors%input_vec, vectors%result_vec, 1.0_dp, &
454 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec)
455 CALL dbcsr_copy(zvec, vectors%result_vec)
456 rsnew = vec_dot_vec(rvec, zvec, mpgrp)
457 beta = rsnew/rsold
458 CALL dbcsr_add(pvec, zvec, beta, 1.0_dp)
459 rsold = rsnew
460 END DO
461 IF (sqrt(rsnew) < control%threshold) control%converged = .true.
462 CALL dbcsr_release(zvec)
463 CALL dbcsr_release(pvec)
464 CALL dbcsr_release(rvec)
465 CALL dbcsr_release(apvec)
466
467 ELSE
468 CALL dbcsr_copy(pvec, xvec)
469 CALL dbcsr_copy(rvec, xvec)
470 CALL dbcsr_set(xvec, 0.0_dp)
471 ! Conjugate gradients
472 rsold = vec_dot_vec(rvec, rvec, mpgrp)
473 DO iter = 1, control%max_iter
474 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_dp, vectors%rep_row_vec, vectors%rep_col_vec)
477 CALL dbcsr_copy(apvec, vectors%result_vec)
478 pap = vec_dot_vec(pvec, apvec, mpgrp)
479 IF (abs(pap) < 1.e-24_dp) THEN
480 alpha = 0.0_dp
481 ELSE
482 alpha = rsold/pap
483 END IF
484 CALL dbcsr_add(xvec, pvec, 1.0_dp, alpha)
485 CALL dbcsr_add(rvec, apvec, 1.0_dp, -alpha)
486 rsnew = vec_dot_vec(rvec, rvec, mpgrp)
487 IF (sqrt(rsnew) < control%threshold) EXIT
488 cpassert(alpha /= 0.0_dp)
489 beta = rsnew/rsold
490 CALL dbcsr_add(pvec, rvec, beta, 1.0_dp)
491 rsold = rsnew
492 END DO
493 IF (sqrt(rsnew) < control%threshold) control%converged = .true.
494 CALL dbcsr_release(pvec)
495 CALL dbcsr_release(rvec)
496 CALL dbcsr_release(apvec)
497 END IF
498
499 END SUBROUTINE conjugate_gradient
500
501! **************************************************************************************************
502!> \brief ...
503!> \param avec ...
504!> \param bvec ...
505!> \param mpgrp ...
506!> \return ...
507! **************************************************************************************************
508 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: av, bv
516 TYPE(dbcsr_iterator_type) :: dbcsr_iter
517
518 adotb = 0.0_dp
519 CALL dbcsr_iterator_start(dbcsr_iter, avec)
520 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
521 CALL dbcsr_iterator_next_block(dbcsr_iter, i, j, av)
522 CALL dbcsr_get_block_p(bvec, i, j, bv, found)
523 IF (found .AND. SIZE(bv) > 0) THEN
524 adotb = adotb + dot_product(av(:, 1), bv(:, 1))
525 END IF
526 END DO
527 CALL dbcsr_iterator_stop(dbcsr_iter)
528 CALL mpgrp%sum(adotb)
529
530 END FUNCTION vec_dot_vec
531
532END MODULE arnoldi_api
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:16
arnoldi_api::arnoldi_ev
subroutine, public arnoldi_ev(matrix, arnoldi_env)
Driver routine for different arnoldi eigenvalue methods the selection which one is to be taken is mad...
Definition arnoldi_api.F:69
arnoldi_api::arnoldi_extremal
subroutine, public arnoldi_extremal(matrix_a, max_ev, min_ev, converged, threshold, max_iter)
simple wrapper to estimate extremal eigenvalues with arnoldi, using the old lanczos interface this hi...
Definition arnoldi_api.F:254
arnoldi_api::arnoldi_conjugate_gradient
subroutine, public arnoldi_conjugate_gradient(matrix_a, vec_x, matrix_p, converged, threshold, max_iter)
Wrapper for conjugated gradient algorithm for Ax=b.
Definition arnoldi_api.F:303
arnoldi_data_methods
The methods which allow to analyze and manipulate the arnoldi procedure The main routine and this sho...
Definition arnoldi_data_methods.F:17
arnoldi_data_methods::deallocate_arnoldi_env
subroutine, public deallocate_arnoldi_env(arnoldi_env)
Deallocate the data in arnoldi_env.
Definition arnoldi_data_methods.F:233
arnoldi_data_methods::get_nrestart
integer function, public get_nrestart(arnoldi_env)
returns the number of restarts allowed for arnoldi
Definition arnoldi_data_methods.F:316
arnoldi_data_methods::arnoldi_is_converged
logical function, public arnoldi_is_converged(arnoldi_env)
Find out whether the method with the current search criterion is converged.
Definition arnoldi_data_methods.F:364
arnoldi_data_methods::get_selected_ritz_vector
subroutine, public get_selected_ritz_vector(arnoldi_env, ind, matrix, vector)
...
Definition arnoldi_data_methods.F:192
arnoldi_data_methods::select_evals
subroutine, public select_evals(arnoldi_env)
perform the selection of eigenvalues, fills the selected_ind array
Definition arnoldi_data_methods.F:258
arnoldi_data_methods::setup_arnoldi_env
subroutine, public setup_arnoldi_env(arnoldi_env, matrix, max_iter, threshold, selection_crit, nval_request, nrestarts, generalized_ev, iram)
This routine sets the environment for the arnoldi iteration and the krylov subspace creation....
Definition arnoldi_data_methods.F:60
arnoldi_data_methods::set_arnoldi_initial_vector
subroutine, public set_arnoldi_initial_vector(arnoldi_env, vector)
...
Definition arnoldi_data_methods.F:433
arnoldi_data_methods::get_selected_ritz_val
complex(dp) function, public get_selected_ritz_val(arnoldi_env, ind)
get a single specific Ritz value from the set of selected
Definition arnoldi_data_methods.F:381
arnoldi_methods
methods for arnoldi iteration
Definition arnoldi_methods.F:16
arnoldi_methods::compute_evals
subroutine, public compute_evals(arnoldi_env)
Call the correct eigensolver, in the arnoldi method only the right eigenvectors are used....
Definition arnoldi_methods.F:153
arnoldi_methods::arnoldi_init
subroutine, public arnoldi_init(matrix, vectors, arnoldi_env)
Interface for the initialization of the arnoldi subspace creation currently it can only setup a rando...
Definition arnoldi_methods.F:198
arnoldi_methods::arnoldi_iram
subroutine, public arnoldi_iram(arnoldi_env)
Alogorithm for the implicit restarts in the arnoldi method this is an early implementation which scal...
Definition arnoldi_methods.F:56
arnoldi_methods::gev_build_subspace
subroutine, public gev_build_subspace(matrix, vectors, arnoldi_env)
builds the basis rothogonal wrt. the metric. The structure looks similar to normal arnoldi but norms,...
Definition arnoldi_methods.F:475
arnoldi_methods::gev_update_data
subroutine, public gev_update_data(matrix, matrix_arnoldi, vectors, arnoldi_env)
Updates all data after an inner loop of the generalized ev arnoldi. Updates rho and C=A-rho*B accordi...
Definition arnoldi_methods.F:585
arnoldi_methods::gev_arnoldi_init
subroutine, public gev_arnoldi_init(matrix, matrix_arnoldi, vectors, arnoldi_env)
Computes the initial guess for the solution of the generalized eigenvalue using the arnoldi method.
Definition arnoldi_methods.F:289
arnoldi_methods::build_subspace
subroutine, public build_subspace(matrix, vectors, arnoldi_env)
Here we create the Krylov subspace and fill the Hessenberg matrix convergence check is only performed...
Definition arnoldi_methods.F:385
arnoldi_types
collection of types used in arnoldi
Definition arnoldi_types.F:16
arnoldi_types::get_control
type(arnoldi_control_type) function, pointer, public get_control(arnoldi_env)
...
Definition arnoldi_types.F:132
arnoldi_vector
operations for skinny matrices/vectors expressed in dbcsr form
Definition arnoldi_vector.F:16
arnoldi_vector::create_replicated_col_vec_from_matrix
subroutine, public create_replicated_col_vec_from_matrix(dbcsr_vec, matrix, ncol)
creates a col vector like object whose blocks can be replicated along the processor row and has the s...
Definition arnoldi_vector.F:163
arnoldi_vector::create_col_vec_from_matrix
subroutine, public create_col_vec_from_matrix(dbcsr_vec, matrix, ncol)
creates a dbcsr col vector like object which lives on proc_col 0 and has the same row dist as the tem...
Definition arnoldi_vector.F:83
arnoldi_vector::create_replicated_row_vec_from_matrix
subroutine, public create_replicated_row_vec_from_matrix(dbcsr_vec, matrix, nrow)
creates a row vector like object whose blocks can be replicated along the processor col and has the s...
Definition arnoldi_vector.F:205
arnoldi_vector::dbcsr_matrix_colvec_multiply
subroutine, public dbcsr_matrix_colvec_multiply(matrix, vec_in, vec_out, alpha, beta, work_row, work_col)
the real driver routine for the multiply, not all symmetries implemented yet
Definition arnoldi_vector.F:531
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
Definition cp_dbcsr_api.F:956
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:930
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1027
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:679
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:794
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:989
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
arnoldi_types::arnoldi_control_type
Definition arnoldi_types.F:26
arnoldi_types::arnoldi_env_type
Definition arnoldi_types.F:54
arnoldi_types::m_x_v_vectors_type
Definition arnoldi_types.F:59
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:186
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
message_passing::mp_comm_type
Definition message_passing.F:189