(git:374b731)
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-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
15MODULE 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
58CONTAINS
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 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 control => get_control(arnoldi_data)
75
76 IF (control%generalized_ev) THEN
77 CALL arnoldi_generalized_ev(matrix, arnoldi_data)
78 ELSE
79 CALL arnoldi_normal_ev(matrix, arnoldi_data)
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_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 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 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_data)
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_data)
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_data, restart_vec)
129 CALL arnoldi_init(matrix, vectors, arnoldi_data)
130 ELSE
131! perform an implicit restart
132 CALL arnoldi_iram(arnoldi_data)
133 END IF
134
135! Generate the subspace
136 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 CALL compute_evals(arnoldi_data)
140
141! Select the evals according to user selection and keep them in arnoldi_data
142 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 IF (.NOT. ASSOCIATED(restart_vec)) ALLOCATE (restart_vec)
146 CALL get_selected_ritz_vector(arnoldi_data, 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_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 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 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_data)
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_data)
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_data)
212 END IF
213
214! Generate the subspace
215 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 CALL compute_evals(arnoldi_data)
219
220! Select the evals according to user selection and keep them in arnoldi_data
221 CALL select_evals(arnoldi_data)
222
223! update the matrices and compute the convergence
224 CALL gev_update_data(matrix, matrix_arnoldi, vectors, arnoldi_data)
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_data_type) :: my_arnoldi
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_data(my_arnoldi, 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, my_arnoldi)
283 converged = arnoldi_is_converged(my_arnoldi)
284 max_ev = real(get_selected_ritz_val(my_arnoldi, 2), dp)
285 min_ev = real(get_selected_ritz_val(my_arnoldi, 1), dp)
286 CALL deallocate_arnoldi_data(my_arnoldi)
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_data_type) :: my_arnoldi
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_data(my_arnoldi, 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(my_arnoldi, arnoldi_matrices, vectors)
358
359 converged = arnoldi_is_converged(my_arnoldi)
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(my_arnoldi)
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_data(my_arnoldi)
380 DEALLOCATE (arnoldi_matrices)
381
382 CALL timestop(handle)
383
384 END SUBROUTINE arnoldi_conjugate_gradient
385
386! **************************************************************************************************
387!> \brief ...
388!> \param arnoldi_data ...
389!> \param arnoldi_matrices ...
390!> \param vectors ...
391! **************************************************************************************************
392 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 control => get_control(arnoldi_data)
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
dbcsr_vector::dbcsr_matrix_colvec_multiply
Definition dbcsr_vector.F:92
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:15
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_api::arnoldi_ev
subroutine, public arnoldi_ev(matrix, arnoldi_data)
Driver routine for different arnoldi eigenvalue methods the selection which one is to be taken is mad...
Definition arnoldi_api.F:69
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::select_evals
subroutine, public select_evals(arnoldi_data)
perform the selection of eigenvalues, fills the selected_ind array
Definition arnoldi_data_methods.F:205
arnoldi_data_methods::setup_arnoldi_data
subroutine, public setup_arnoldi_data(arnoldi_data, 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:64
arnoldi_data_methods::get_selected_ritz_val
complex(real_8) function, public get_selected_ritz_val(ar_data, ind)
get a single specific Ritz value from the set of selected
Definition arnoldi_data_methods.F:299
arnoldi_data_methods::get_selected_ritz_vector
subroutine, public get_selected_ritz_vector(arnoldi_data, ind, matrix, vector)
Deallocate the data in arnoldi_data.
Definition arnoldi_data_methods.F:173
arnoldi_data_methods::deallocate_arnoldi_data
subroutine, public deallocate_arnoldi_data(arnoldi_data)
Deallocate the data in arnoldi_data.
Definition arnoldi_data_methods.F:187
arnoldi_data_methods::arnoldi_is_converged
logical function, public arnoldi_is_converged(ar_data)
Find out whether the method with the current search criterion is converged.
Definition arnoldi_data_methods.F:282
arnoldi_data_methods::set_arnoldi_initial_vector
subroutine, public set_arnoldi_initial_vector(ar_data, vector)
...
Definition arnoldi_data_methods.F:355
arnoldi_data_methods::get_nrestart
integer function, public get_nrestart(arnoldi_data)
returns the number of restarts allowed for arnoldi
Definition arnoldi_data_methods.F:233
arnoldi_methods
methods for arnoldi iteration
Definition arnoldi_methods.F:15
arnoldi_methods::arnoldi_init
subroutine, public arnoldi_init(matrix, vectors, arnoldi_data)
Interface for the initialization of the arnoldi subspace creation currently it can only setup a rando...
Definition arnoldi_methods.F:100
arnoldi_methods::gev_update_data
subroutine, public gev_update_data(matrix, matrix_arnoldi, vectors, arnoldi_data)
in the generalized eigenvalue the matrix depende on the projection therefore the outer loop has to bu...
Definition arnoldi_methods.F:202
arnoldi_methods::gev_build_subspace
subroutine, public gev_build_subspace(matrix, vectors, arnoldi_data)
here the iterations are performed and the krylov space for the generalized eigenvalue problem is crea...
Definition arnoldi_methods.F:175
arnoldi_methods::compute_evals
subroutine, public compute_evals(arnoldi_data)
Interface to compute the eigenvalues of a nonsymmetric matrix This is only the serial version.
Definition arnoldi_methods.F:76
arnoldi_methods::gev_arnoldi_init
subroutine, public gev_arnoldi_init(matrix, matrix_arnoldi, vectors, arnoldi_data)
Interface for the initialization of the arnoldi subspace creation for the generalized eigenvalue prob...
Definition arnoldi_methods.F:126
arnoldi_methods::build_subspace
subroutine, public build_subspace(matrix, vectors, arnoldi_data)
here the iterations are performed and the krylov space is constructed
Definition arnoldi_methods.F:150
arnoldi_methods::arnoldi_iram
subroutine, public arnoldi_iram(arnoldi_data)
Interface for the routine calcualting the implicit restarts Currently all based on lapack.
Definition arnoldi_methods.F:55
arnoldi_types
collection of types used in arnoldi
Definition arnoldi_types.F:15
arnoldi_types::get_control
type(arnoldi_control_type) function, pointer, public get_control(ar_data)
...
Definition arnoldi_types.F:245
dbcsr_vector
operations for skinny matrices/vectors expressed in dbcsr form
Definition dbcsr_vector.F:15
dbcsr_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 dbcsr_vector.F:236
dbcsr_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 dbcsr_vector.F:109
dbcsr_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 dbcsr_vector.F:193
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_data_type
Definition arnoldi_types.F:84
arnoldi_types::m_x_v_vectors_type
Definition arnoldi_types.F:92
message_passing::mp_comm_type
Definition message_passing.F:189