(git:374b731)
Loading...
Searching...
No Matches
arnoldi_methods.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 methods for arnoldi iteration
10!> \par History
11!> 2014.09 created [Florian Schiffmann]
12!> \author Florian Schiffmann
13! **************************************************************************************************
14
15MODULE arnoldi_methods
16 USE arnoldi_geev, ONLY: arnoldi_general_local_diag, &
17 arnoldi_symm_local_diag, &
18 arnoldi_tridiag_local_diag
19 USE arnoldi_types, ONLY: &
20 arnoldi_control_type, arnoldi_data_c_type, arnoldi_data_d_type, arnoldi_data_s_type, &
21 arnoldi_data_type, arnoldi_data_z_type, get_control, get_data_c, get_data_d, get_data_s, &
22 get_data_z, has_d_cmplx, has_d_real, m_x_v_vectors_type
23 USE dbcsr_api, ONLY: &
24 dbcsr_add, dbcsr_copy, dbcsr_get_data_p, dbcsr_get_info, dbcsr_iterator_blocks_left, &
25 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
26 dbcsr_p_type, dbcsr_scale, dbcsr_type
27 USE dbcsr_vector, ONLY: dbcsr_matrix_colvec_multiply
28 USE kinds, ONLY: real_8
29 USE message_passing, ONLY: mp_comm_type
30#include "../base/base_uses.f90"
31
32 IMPLICIT NONE
33
34 PRIVATE
35
36 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'arnoldi_methods'
37
38 PUBLIC :: arnoldi_init, build_subspace, compute_evals, arnoldi_iram, &
39 gev_arnoldi_init, gev_build_subspace, gev_update_data
40
41 INTERFACE convert_matrix
42 MODULE PROCEDURE convert_matrix_z_to_d
43 MODULE PROCEDURE convert_matrix_d_to_z
44 MODULE PROCEDURE convert_matrix_z_to_z
45 END INTERFACE
46
47CONTAINS
48
49! **************************************************************************************************
50!> \brief Interface for the routine calcualting the implicit restarts
51!> Currently all based on lapack
52!> \param arnoldi_data ...
53! **************************************************************************************************
54 SUBROUTINE arnoldi_iram(arnoldi_data)
55 TYPE(arnoldi_data_type) :: arnoldi_data
56
57 CHARACTER(LEN=*), PARAMETER :: routinen = 'arnoldi_iram'
58
59 INTEGER :: handle
60
61 CALL timeset(routinen, handle)
62
63 IF (has_d_real(arnoldi_data)) CALL arnoldi_iram_d(arnoldi_data)
64 IF (has_d_cmplx(arnoldi_data)) CALL arnoldi_iram_z(arnoldi_data)
65
66 CALL timestop(handle)
67
68 END SUBROUTINE arnoldi_iram
69
70! **************************************************************************************************
71!> \brief Interface to compute the eigenvalues of a nonsymmetric matrix
72!> This is only the serial version
73!> \param arnoldi_data ...
74! **************************************************************************************************
75 SUBROUTINE compute_evals(arnoldi_data)
76 TYPE(arnoldi_data_type) :: arnoldi_data
77
78 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_evals'
79
80 INTEGER :: handle
81
82 CALL timeset(routinen, handle)
83
84 IF (has_d_real(arnoldi_data)) CALL compute_evals_d(arnoldi_data)
85 IF (has_d_cmplx(arnoldi_data)) CALL compute_evals_z(arnoldi_data)
86
87 CALL timestop(handle)
88
89 END SUBROUTINE compute_evals
90
91! **************************************************************************************************
92!> \brief Interface for the initialization of the arnoldi subspace creation
93!> currently it can only setup a random vector but can be improved to
94!> various types of restarts easily
95!> \param matrix pointer to the matrices as described in main interface
96!> \param vectors work vectors for the matrix vector multiplications
97!> \param arnoldi_data all data concerning the subspace
98! **************************************************************************************************
99 SUBROUTINE arnoldi_init(matrix, vectors, arnoldi_data)
100 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
101 TYPE(m_x_v_vectors_type) :: vectors
102 TYPE(arnoldi_data_type) :: arnoldi_data
103
104 CHARACTER(LEN=*), PARAMETER :: routinen = 'arnoldi_init'
105
106 INTEGER :: handle
107
108 CALL timeset(routinen, handle)
109
110 IF (has_d_real(arnoldi_data)) CALL arnoldi_init_d(matrix, vectors, arnoldi_data)
111 IF (has_d_cmplx(arnoldi_data)) CALL arnoldi_init_z(matrix, vectors, arnoldi_data)
112
113 CALL timestop(handle)
114
115 END SUBROUTINE arnoldi_init
116
117! **************************************************************************************************
118!> \brief Interface for the initialization of the arnoldi subspace creation
119!> for the generalized eigenvalue problem
120!> \param matrix pointer to the matrices as described in main interface
121!> \param matrix_arnoldi ...
122!> \param vectors work vectors for the matrix vector multiplications
123!> \param arnoldi_data all data concerning the subspace
124! **************************************************************************************************
125 SUBROUTINE gev_arnoldi_init(matrix, matrix_arnoldi, vectors, arnoldi_data)
126 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix, matrix_arnoldi
127 TYPE(m_x_v_vectors_type) :: vectors
128 TYPE(arnoldi_data_type) :: arnoldi_data
129
130 CHARACTER(LEN=*), PARAMETER :: routinen = 'gev_arnoldi_init'
131
132 INTEGER :: handle
133
134 CALL timeset(routinen, handle)
135
136 IF (has_d_real(arnoldi_data)) CALL gev_arnoldi_init_d(matrix, matrix_arnoldi, vectors, arnoldi_data)
137 IF (has_d_cmplx(arnoldi_data)) CALL gev_arnoldi_init_z(matrix, matrix_arnoldi, vectors, arnoldi_data)
138
139 CALL timestop(handle)
140
141 END SUBROUTINE gev_arnoldi_init
142
143! **************************************************************************************************
144!> \brief here the iterations are performed and the krylov space is constructed
145!> \param matrix see above
146!> \param vectors see above
147!> \param arnoldi_data see above
148! **************************************************************************************************
149 SUBROUTINE build_subspace(matrix, vectors, arnoldi_data)
150 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
151 TYPE(m_x_v_vectors_type) :: vectors
152 TYPE(arnoldi_data_type) :: arnoldi_data
153
154 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_subspace'
155
156 INTEGER :: handle
157
158 CALL timeset(routinen, handle)
159
160 IF (has_d_real(arnoldi_data)) CALL build_subspace_d(matrix, vectors, arnoldi_data)
161 IF (has_d_cmplx(arnoldi_data)) CALL build_subspace_z(matrix, vectors, arnoldi_data)
162
163 CALL timestop(handle)
164
165 END SUBROUTINE build_subspace
166
167! **************************************************************************************************
168!> \brief here the iterations are performed and the krylov space for the generalized
169!> eigenvalue problem is created
170!> \param matrix see above
171!> \param vectors see above
172!> \param arnoldi_data see above
173! **************************************************************************************************
174 SUBROUTINE gev_build_subspace(matrix, vectors, arnoldi_data)
175 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
176 TYPE(m_x_v_vectors_type) :: vectors
177 TYPE(arnoldi_data_type) :: arnoldi_data
178
179 CHARACTER(LEN=*), PARAMETER :: routinen = 'gev_build_subspace'
180
181 INTEGER :: handle
182
183 CALL timeset(routinen, handle)
184
185 IF (has_d_real(arnoldi_data)) CALL gev_build_subspace_d(matrix, vectors, arnoldi_data)
186 IF (has_d_cmplx(arnoldi_data)) CALL gev_build_subspace_z(matrix, vectors, arnoldi_data)
187
188 CALL timestop(handle)
189
190 END SUBROUTINE gev_build_subspace
191
192! **************************************************************************************************
193!> \brief in the generalized eigenvalue the matrix depende on the projection
194!> therefore the outer loop has to build a new set of matrices for the
195!> inner loop
196!> \param matrix see above
197!> \param matrix_arnoldi ...
198!> \param vectors ...
199!> \param arnoldi_data see above
200! **************************************************************************************************
201 SUBROUTINE gev_update_data(matrix, matrix_arnoldi, vectors, arnoldi_data)
202 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix, matrix_arnoldi
203 TYPE(m_x_v_vectors_type) :: vectors
204 TYPE(arnoldi_data_type) :: arnoldi_data
205
206 CHARACTER(LEN=*), PARAMETER :: routinen = 'gev_update_data'
207
208 INTEGER :: handle
209
210 CALL timeset(routinen, handle)
211
212 IF (has_d_real(arnoldi_data)) CALL gev_update_data_d(matrix, matrix_arnoldi, vectors, arnoldi_data)
213 IF (has_d_cmplx(arnoldi_data)) CALL gev_update_data_z(matrix, matrix_arnoldi, vectors, arnoldi_data)
214
215 CALL timestop(handle)
216
217 END SUBROUTINE gev_update_data
218
219! **************************************************************************************************
220!> \brief ...
221!> \param m_out ...
222!> \param m_in ...
223! **************************************************************************************************
224 SUBROUTINE convert_matrix_z_to_d(m_out, m_in)
225 REAL(real_8), DIMENSION(:, :) :: m_out
226 COMPLEX(real_8), DIMENSION(:, :) :: m_in
227
228 m_out(:, :) = real(m_in(:, :), kind=real_8)
229 END SUBROUTINE convert_matrix_z_to_d
230
231! **************************************************************************************************
232!> \brief ...
233!> \param m_out ...
234!> \param m_in ...
235! **************************************************************************************************
236 SUBROUTINE convert_matrix_d_to_z(m_out, m_in)
237 COMPLEX(real_8), DIMENSION(:, :) :: m_out
238 REAL(real_8), DIMENSION(:, :) :: m_in
239
240 m_out(:, :) = cmplx(m_in, 0.0, kind=real_8)
241 END SUBROUTINE convert_matrix_d_to_z
242
243! I kno that one is stupid but like this it simplifies the template
244! **************************************************************************************************
245!> \brief ...
246!> \param m_out ...
247!> \param m_in ...
248! **************************************************************************************************
249 SUBROUTINE convert_matrix_z_to_z(m_out, m_in)
250 COMPLEX(real_8), DIMENSION(:, :) :: m_out, m_in
251
252 m_out(:, :) = m_in
253 END SUBROUTINE convert_matrix_z_to_z
254
255! **************************************************************************************************
256!> \brief Call the correct eigensolver, in the arnoldi method only the right
257!> eigenvectors are used. Lefts are created here but dumped immediately
258!> \param arnoldi_data ...
259! **************************************************************************************************
260 SUBROUTINE compute_evals_d (arnoldi_data)
261 TYPE(arnoldi_data_type) :: arnoldi_data
262
263 COMPLEX(real_8), DIMENSION(:, :), ALLOCATABLE :: levec
264 TYPE(arnoldi_data_d_type), POINTER :: ar_data
265 INTEGER :: ndim
266 TYPE(arnoldi_control_type), POINTER :: control
267
268 ar_data => get_data_d(arnoldi_data)
269 control => get_control(arnoldi_data)
270 ndim = control%current_step
271 ALLOCATE (levec(ndim, ndim))
272
273! Needs antoher interface as the calls to real and complex geev differ (sucks!)
274! only perform the diagonalization on processors which hold data
275 IF (control%generalized_ev) THEN
276 CALL arnoldi_symm_local_diag('V', ar_data%Hessenberg(1:ndim, 1:ndim), ndim, &
277 ar_data%evals(1:ndim), ar_data%revec(1:ndim, 1:ndim))
278 ELSE
279 IF (control%symmetric) THEN
280 CALL arnoldi_tridiag_local_diag('N', 'V', ar_data%Hessenberg(1:ndim, 1:ndim), ndim, &
281 ar_data%evals(1:ndim), ar_data%revec(1:ndim, 1:ndim), levec)
282 ELSE
283 CALL arnoldi_general_local_diag('N', 'V', ar_data%Hessenberg(1:ndim, 1:ndim), ndim, &
284 ar_data%evals(1:ndim), ar_data%revec(1:ndim, 1:ndim), levec)
285 END IF
286 END IF
287
288 DEALLOCATE (levec)
289
290 END SUBROUTINE compute_evals_d
291
292! **************************************************************************************************
293!> \brief Initialization of the arnoldi vector. Here a random vector is used,
294!> might be generalized in the future
295!> \param matrix ...
296!> \param vectors ...
297!> \param arnoldi_data ...
298! **************************************************************************************************
299 SUBROUTINE arnoldi_init_d (matrix, vectors, arnoldi_data)
300 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
301 TYPE(m_x_v_vectors_type) :: vectors
302 TYPE(arnoldi_data_type) :: arnoldi_data
303
304 INTEGER :: i, iseed(4), row_size, col_size, &
305 nrow_local, ncol_local, &
306 row, col
307 REAL(real_8) :: rnorm
308 TYPE(dbcsr_iterator_type) :: iter
309 REAL(kind=real_8) :: norm
310 REAL(kind=real_8), DIMENSION(:), ALLOCATABLE :: &
311 v_vec, w_vec
312 REAL(kind=real_8), DIMENSION(:), POINTER :: data_vec
313 LOGICAL :: local_comp
314 TYPE(arnoldi_data_d_type), POINTER :: ar_data
315 TYPE(arnoldi_control_type), POINTER :: control
316 TYPE(mp_comm_type) :: pcol_group
317
318 control => get_control(arnoldi_data)
319 pcol_group = control%pcol_group
320 local_comp = control%local_comp
321
322 ar_data => get_data_d(arnoldi_data)
323
324 ! create a local data copy to store the vectors and make Gram Schmidt a bit simpler
325 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
326 ALLOCATE (v_vec(nrow_local))
327 ALLOCATE (w_vec(nrow_local))
328 v_vec = 0.0_real_8; w_vec = 0.0_real_8
329 ar_data%Hessenberg = 0.0_real_8
330
331 IF (control%has_initial_vector) THEN
332 ! after calling the set routine the initial vector is stored in f_vec
333 CALL transfer_local_array_to_dbcsr_d (vectors%input_vec, ar_data%f_vec, nrow_local, control%local_comp)
334 ELSE
335 ! Setup the initial normalized random vector (sufficient if it only happens on proc_col 0)
336 CALL dbcsr_iterator_start(iter, vectors%input_vec)
337 DO WHILE (dbcsr_iterator_blocks_left(iter))
338 CALL dbcsr_iterator_next_block(iter, row, col, data_vec, row_size=row_size, col_size=col_size)
339 iseed(1) = 2; iseed(2) = mod(row, 4095); iseed(3) = mod(col, 4095); iseed(4) = 11
340 CALL dlarnv(2, iseed, row_size*col_size, data_vec)
341 END DO
342 CALL dbcsr_iterator_stop(iter)
343 END IF
344
345 CALL transfer_dbcsr_to_local_array_d (vectors%input_vec, v_vec, nrow_local, control%local_comp)
346
347 ! compute the vector norm of the random vectorm, get it real valued as well (rnorm)
348 CALL compute_norms_d (v_vec, norm, rnorm, control%pcol_group)
349
350 IF (rnorm == 0) rnorm = 1 ! catch case where this rank has no actual data
351 CALL dbcsr_scale(vectors%input_vec, real(1.0, real_8)/rnorm)
352
353 ! Everything prepared, initialize the Arnoldi iteration
354 CALL transfer_dbcsr_to_local_array_d (vectors%input_vec, v_vec, nrow_local, control%local_comp)
355
356 ! This permits to compute the subspace of a matrix which is a product of multiple matrices
357 DO i = 1, SIZE(matrix)
358 CALL dbcsr_matrix_colvec_multiply(matrix(i)%matrix, vectors%input_vec, vectors%result_vec, 1.0_real_8, &
359 0.0_real_8, vectors%rep_row_vec, vectors%rep_col_vec)
360 CALL dbcsr_copy(vectors%input_vec, vectors%result_vec)
361 END DO
362
363 CALL transfer_dbcsr_to_local_array_d (vectors%result_vec, w_vec, nrow_local, control%local_comp)
364
365 ! Put the projection into the Hessenberg matrix, and make the vectors orthonormal
366 ar_data%Hessenberg(1, 1) = dot_product(v_vec, w_vec)
367 CALL pcol_group%sum(ar_data%Hessenberg(1, 1))
368 ar_data%f_vec = w_vec - v_vec*ar_data%Hessenberg(1, 1)
369
370 ar_data%local_history(:, 1) = v_vec(:)
371
372 ! We did the first step in here so we should set the current step for the subspace generation accordingly
373 control%current_step = 1
374
375 DEALLOCATE (v_vec, w_vec)
376
377 END SUBROUTINE arnoldi_init_d
378
379! **************************************************************************************************
380!> \brief Alogorithm for the implicit restarts in the arnoldi method
381!> this is an early implementation which scales subspace size^4
382!> by replacing the lapack calls with direct math the
383!> QR and gemms can be made linear and a N^2 sacling will be acchieved
384!> however this already sets the framework but should be used with care
385!> \param arnoldi_data ...
386! **************************************************************************************************
387 SUBROUTINE arnoldi_iram_d (arnoldi_data)
388 TYPE(arnoldi_data_type) :: arnoldi_data
389
390 TYPE(arnoldi_data_d_type), POINTER :: ar_data
391 TYPE(arnoldi_control_type), POINTER :: control
392 COMPLEX(real_8), DIMENSION(:, :), ALLOCATABLE :: tmp_mat, safe_mat, q, tmp_mat1
393 COMPLEX(real_8), DIMENSION(:), ALLOCATABLE :: work, tau, work_measure
394 INTEGER :: msize, lwork, i, j, info, nwant
395 REAL(kind=real_8) :: beta, sigma
396 REAL(kind=real_8), DIMENSION(:, :), ALLOCATABLE :: qdata
397
398 ! This is just a terribly inefficient implementation but I hope it is correct and might serve as a reference
399 ar_data => get_data_d(arnoldi_data)
400 control => get_control(arnoldi_data)
401 msize = control%current_step
402 nwant = control%nval_out
403 ALLOCATE (tmp_mat(msize, msize)); ALLOCATE (safe_mat(msize, msize))
404 ALLOCATE (q(msize, msize)); ALLOCATE (tmp_mat1(msize, msize))
405 ALLOCATE (work_measure(1))
406 ALLOCATE (tau(msize)); ALLOCATE (qdata(msize, msize))
407 !make Q identity
408 q = cmplx(0.0, 0.0, real_8)
409 DO i = 1, msize
410 q(i, i) = cmplx(1.0, 0.0, real_8)
411 END DO
412
413 ! Looks a bit odd, but safe_mat will contain the result in the end, while tmpmat gets violated by lapack
414 CALL convert_matrix(tmp_mat, ar_data%Hessenberg(1:msize, 1:msize))
415 safe_mat(:, :) = tmp_mat(:, :)
416
417 DO i = 1, msize
418 ! A bit a strange check but in the end we only want to shift the unwanted evals
419 IF (any(control%selected_ind == i)) cycle
420 ! Here we shift the matrix by subtracting unwanted evals from the diagonal
421 DO j = 1, msize
422 tmp_mat(j, j) = tmp_mat(j, j) - ar_data%evals(i)
423 END DO
424 ! Now we repair the damage by QR factorizing
425 lwork = -1
426 CALL zgeqrf(msize, msize, tmp_mat, msize, tau, work_measure, lwork, info)
427 lwork = int(work_measure(1))
428 IF (ALLOCATED(work)) THEN
429 IF (SIZE(work) .LT. lwork) THEN
430 DEALLOCATE (work)
431 END IF
432 END IF
433 IF (.NOT. ALLOCATED(work)) ALLOCATE (work(lwork))
434 CALL zgeqrf(msize, msize, tmp_mat, msize, tau, work, lwork, info)
435 ! Ask Lapack to reconstruct Q from its own way of storing data (tmpmat will contain Q)
436 CALL zungqr(msize, msize, msize, tmp_mat, msize, tau, work, lwork, info)
437 ! update Q=Q*Q_current
438 tmp_mat1(:, :) = q(:, :)
439 CALL zgemm('N', 'N', msize, msize, msize, cmplx(1.0, 0.0, real_8), tmp_mat1, msize, tmp_mat, msize, cmplx(0.0, 0.0,&
440 & real_8), &
441 q, msize)
442 ! Update H=(Q*)HQ
443 CALL zgemm('C', 'N', msize, msize, msize, cmplx(1.0, 0.0, real_8), tmp_mat, msize, safe_mat, msize, cmplx(0.0, 0.0,&
444 & real_8), &
445 tmp_mat1, msize)
446 CALL zgemm('N', 'N', msize, msize, msize, cmplx(1.0, 0.0, real_8), tmp_mat1, msize, tmp_mat, msize, cmplx(0.0, 0.0,&
447 & real_8), &
448 safe_mat, msize)
449
450 ! this one is crucial for numerics not to accumulate noise in the subdiagonals
451 DO j = 1, msize
452 safe_mat(j + 2:msize, j) = cmplx(0.0, 0.0, real_8)
453 END DO
454 tmp_mat(:, :) = safe_mat(:, :)
455 END DO
456
457 ! Now we can compute our restart quantities
458 ar_data%Hessenberg = 0.0_real_8
459 CALL convert_matrix(ar_data%Hessenberg(1:msize, 1:msize), safe_mat)
460 CALL convert_matrix(qdata, q)
461
462 beta = ar_data%Hessenberg(nwant + 1, nwant); sigma = qdata(msize, nwant)
463
464 !update the residuum and the basis vectors
465 IF (control%local_comp) THEN
466 ar_data%f_vec = matmul(ar_data%local_history(:, 1:msize), qdata(1:msize, nwant + 1))*beta + ar_data%f_vec(:)*sigma
467 ar_data%local_history(:, 1:nwant) = matmul(ar_data%local_history(:, 1:msize), qdata(1:msize, 1:nwant))
468 END IF
469 ! Set the current step to nwant so the subspace build knows where to start
470 control%current_step = nwant
471
472 DEALLOCATE (tmp_mat, safe_mat, q, qdata, tmp_mat1, work, tau, work_measure)
473
474 END SUBROUTINE arnoldi_iram_d
475
476! **************************************************************************************************
477!> \brief Here we create the Krylov subspace and fill the Hessenberg matrix
478!> convergence check is only performed on subspace convergence
479!> Gram Schidt is used to orthonogonalize.
480!> If this is numericall not sufficient a Daniel, Gragg, Kaufman and Steward
481!> correction is performed
482!> \param matrix ...
483!> \param vectors ...
484!> \param arnoldi_data ...
485! **************************************************************************************************
486 SUBROUTINE build_subspace_d (matrix, vectors, arnoldi_data)
487 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
488 TYPE(m_x_v_vectors_type), TARGET :: vectors
489 TYPE(arnoldi_data_type) :: arnoldi_data
490
491 INTEGER :: i, j, ncol_local, nrow_local
492 REAL(real_8) :: rnorm
493 TYPE(arnoldi_control_type), POINTER :: control
494 TYPE(arnoldi_data_d_type), POINTER :: ar_data
495 REAL(kind=real_8) :: norm
496 REAL(kind=real_8), ALLOCATABLE, DIMENSION(:) :: h_vec, s_vec, v_vec, w_vec
497 TYPE(dbcsr_type), POINTER :: input_vec, result_vec, swap_vec
498
499 ar_data => get_data_d(arnoldi_data)
500 control => get_control(arnoldi_data)
501 control%converged = .false.
502
503 ! create the vectors required during the iterations
504 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
505 ALLOCATE (v_vec(nrow_local)); ALLOCATE (w_vec(nrow_local))
506 v_vec = 0.0_real_8; w_vec = 0.0_real_8
507 ALLOCATE (s_vec(control%max_iter)); ALLOCATE (h_vec(control%max_iter))
508
509 DO j = control%current_step, control%max_iter - 1
510
511 ! compute the vector norm of the residuum, get it real valued as well (rnorm)
512 CALL compute_norms_d (ar_data%f_vec, norm, rnorm, control%pcol_group)
513
514 ! check convergence and inform everybody about it, a bit annoying to talk to everybody because of that
515 IF (control%myproc == 0) control%converged = rnorm .LT. real(control%threshold, real_8)
516 CALL control%mp_group%bcast(control%converged, 0)
517 IF (control%converged) EXIT
518
519 ! transfer normalized residdum to history and its norm to the Hessenberg matrix
520 IF (rnorm == 0) rnorm = 1 ! catch case where this rank has no actual data
521 v_vec(:) = ar_data%f_vec(:)/rnorm; ar_data%local_history(:, j + 1) = v_vec(:); ar_data%Hessenberg(j + 1, j) = norm
522
523 input_vec => vectors%input_vec
524 result_vec => vectors%result_vec
525 CALL transfer_local_array_to_dbcsr_d (input_vec, v_vec, nrow_local, control%local_comp)
526
527 ! This permits to compute the subspace of a matrix which is a product of two matrices
528 DO i = 1, SIZE(matrix)
529 CALL dbcsr_matrix_colvec_multiply(matrix(i)%matrix, input_vec, result_vec, 1.0_real_8, &
530 0.0_real_8, vectors%rep_row_vec, vectors%rep_col_vec)
531 swap_vec => input_vec
532 input_vec => result_vec
533 result_vec => swap_vec
534 END DO
535
536 CALL transfer_dbcsr_to_local_array_d (input_vec, w_vec, nrow_local, control%local_comp)
537
538 ! Let's do the orthonormalization, to get the new f_vec. First try the Gram Schmidt scheme
539 CALL gram_schmidt_ortho_d (h_vec, ar_data%f_vec, s_vec, w_vec, nrow_local, j + 1, &
540 ar_data%local_history, ar_data%local_history, control%local_comp, control%pcol_group)
541
542 ! A bit more expensive but simply always top up with a DGKS correction, otherwise numerics
543 ! can become a problem later on, there is probably a good check whether it's necessary, but we don't perform it
544 CALL dgks_ortho_d (h_vec, ar_data%f_vec, s_vec, nrow_local, j + 1, ar_data%local_history, &
545 ar_data%local_history, control%local_comp, control%pcol_group)
546 ! Finally we can put the projections into our Hessenberg matrix
547 ar_data%Hessenberg(1:j + 1, j + 1) = h_vec(1:j + 1)
548 control%current_step = j + 1
549 END DO
550
551 ! compute the vector norm of the final residuum and put it in to Hessenberg
552 CALL compute_norms_d (ar_data%f_vec, norm, rnorm, control%pcol_group)
553 ar_data%Hessenberg(control%current_step + 1, control%current_step) = norm
554
555 ! broadcast the Hessenberg matrix so we don't need to care later on
556 CALL control%mp_group%bcast(ar_data%Hessenberg, 0)
557
558 DEALLOCATE (v_vec, w_vec, h_vec, s_vec)
559
560 END SUBROUTINE build_subspace_d
561
562! **************************************************************************************************
563!> \brief Helper routine to transfer the all data of a dbcsr matrix to a local array
564!> \param vec ...
565!> \param array ...
566!> \param n ...
567!> \param is_local ...
568! **************************************************************************************************
569 SUBROUTINE transfer_dbcsr_to_local_array_d (vec, array, n, is_local)
570 TYPE(dbcsr_type) :: vec
571 REAL(kind=real_8), DIMENSION(:) :: array
572 INTEGER :: n
573 LOGICAL :: is_local
574 REAL(kind=real_8), DIMENSION(:), POINTER :: data_vec
575
576 data_vec => dbcsr_get_data_p(vec, select_data_type=0.0_real_8)
577 IF (is_local) array(1:n) = data_vec(1:n)
578
579 END SUBROUTINE transfer_dbcsr_to_local_array_d
580
581! **************************************************************************************************
582!> \brief The inverse routine transferring data back from an array to a dbcsr
583!> \param vec ...
584!> \param array ...
585!> \param n ...
586!> \param is_local ...
587! **************************************************************************************************
588 SUBROUTINE transfer_local_array_to_dbcsr_d (vec, array, n, is_local)
589 TYPE(dbcsr_type) :: vec
590 REAL(kind=real_8), DIMENSION(:) :: array
591 INTEGER :: n
592 LOGICAL :: is_local
593 REAL(kind=real_8), DIMENSION(:), POINTER :: data_vec
594
595 data_vec => dbcsr_get_data_p(vec, select_data_type=0.0_real_8)
596 IF (is_local) data_vec(1:n) = array(1:n)
597
598 END SUBROUTINE transfer_local_array_to_dbcsr_d
599
600! **************************************************************************************************
601!> \brief Gram-Schmidt in matrix vector form
602!> \param h_vec ...
603!> \param f_vec ...
604!> \param s_vec ...
605!> \param w_vec ...
606!> \param nrow_local ...
607!> \param j ...
608!> \param local_history ...
609!> \param reorth_mat ...
610!> \param local_comp ...
611!> \param pcol_group ...
612! **************************************************************************************************
613 SUBROUTINE gram_schmidt_ortho_d (h_vec, f_vec, s_vec, w_vec, nrow_local, &
614 j, local_history, reorth_mat, local_comp, pcol_group)
615 REAL(kind=real_8), DIMENSION(:) :: h_vec, s_vec, f_vec, w_vec
616 REAL(kind=real_8), DIMENSION(:, :) :: local_history, reorth_mat
617 INTEGER :: nrow_local, j
618 TYPE(mp_comm_type), INTENT(IN) :: pcol_group
619 LOGICAL :: local_comp
620
621 CHARACTER(LEN=*), PARAMETER :: routinen = 'Gram_Schmidt_ortho_d'
622 INTEGER :: handle
623
624 CALL timeset(routinen, handle)
625
626 ! Let's do the orthonormalization, first try the Gram Schmidt scheme
627 h_vec = 0.0_real_8; f_vec = 0.0_real_8; s_vec = 0.0_real_8
628 IF (local_comp) CALL dgemv('T', nrow_local, j, 1.0_real_8, local_history, &
629 nrow_local, w_vec, 1, 0.0_real_8, h_vec, 1)
630 CALL pcol_group%sum(h_vec(1:j))
631
632 IF (local_comp) CALL dgemv('N', nrow_local, j, 1.0_real_8, reorth_mat, &
633 nrow_local, h_vec, 1, 0.0_real_8, f_vec, 1)
634 f_vec(:) = w_vec(:) - f_vec(:)
635
636 CALL timestop(handle)
637
638 END SUBROUTINE gram_schmidt_ortho_d
639
640! **************************************************************************************************
641!> \brief Compute the Daniel, Gragg, Kaufman and Steward correction
642!> \param h_vec ...
643!> \param f_vec ...
644!> \param s_vec ...
645!> \param nrow_local ...
646!> \param j ...
647!> \param local_history ...
648!> \param reorth_mat ...
649!> \param local_comp ...
650!> \param pcol_group ...
651! **************************************************************************************************
652 SUBROUTINE dgks_ortho_d (h_vec, f_vec, s_vec, nrow_local, j, &
653 local_history, reorth_mat, local_comp, pcol_group)
654 REAL(kind=real_8), DIMENSION(:) :: h_vec, s_vec, f_vec
655 REAL(kind=real_8), DIMENSION(:, :) :: local_history, reorth_mat
656 INTEGER :: nrow_local, j
657 TYPE(mp_comm_type), INTENT(IN) :: pcol_group
658
659 CHARACTER(LEN=*), PARAMETER :: routinen = 'DGKS_ortho_d'
660 LOGICAL :: local_comp
661 INTEGER :: handle
662
663 CALL timeset(routinen, handle)
664
665 IF (local_comp) CALL dgemv('T', nrow_local, j, 1.0_real_8, local_history, &
666 nrow_local, f_vec, 1, 0.0_real_8, s_vec, 1)
667 CALL pcol_group%sum(s_vec(1:j))
668 IF (local_comp) CALL dgemv('N', nrow_local, j, -1.0_real_8, reorth_mat, &
669 nrow_local, s_vec, 1, 1.0_real_8, f_vec, 1)
670 h_vec(1:j) = h_vec(1:j) + s_vec(1:j)
671
672 CALL timestop(handle)
673
674 END SUBROUTINE dgks_ortho_d
675
676! **************************************************************************************************
677!> \brief Compute the norm of a vector distributed along proc_col
678!> as local arrays. Always return the real part next to the complex rep.
679!> \param vec ...
680!> \param norm ...
681!> \param rnorm ...
682!> \param pcol_group ...
683! **************************************************************************************************
684 SUBROUTINE compute_norms_d (vec, norm, rnorm, pcol_group)
685 REAL(kind=real_8), DIMENSION(:) :: vec
686 REAL(real_8) :: rnorm
687 REAL(kind=real_8) :: norm
688 TYPE(mp_comm_type), INTENT(IN) :: pcol_group
689
690 ! the norm is computed along the processor column
691 norm = dot_product(vec, vec)
692 CALL pcol_group%sum(norm)
693 rnorm = sqrt(real(norm, real_8))
694 norm=rnorm
695
696 END SUBROUTINE compute_norms_d
697
698! **************************************************************************************************
699!> \brief Computes the initial guess for the solution of the generalized eigenvalue
700!> using the arnoldi method
701!> \param matrix ...
702!> \param matrix_arnoldi ...
703!> \param vectors ...
704!> \param arnoldi_data ...
705! **************************************************************************************************
706
707 SUBROUTINE gev_arnoldi_init_d (matrix, matrix_arnoldi, vectors, arnoldi_data)
708 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
709 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix_arnoldi
710 TYPE(m_x_v_vectors_type) :: vectors
711 TYPE(arnoldi_data_type) :: arnoldi_data
712
713 INTEGER :: iseed(4), row_size, col_size, &
714 nrow_local, ncol_local, &
715 row, col
716 REAL(real_8) :: rnorm
717 TYPE(dbcsr_iterator_type) :: iter
718 REAL(kind=real_8) :: norm, denom
719 REAL(kind=real_8), DIMENSION(:), ALLOCATABLE :: &
720 v_vec, w_vec
721 REAL(kind=real_8), DIMENSION(:), POINTER :: data_vec
722 LOGICAL :: local_comp
723 TYPE(arnoldi_data_d_type), POINTER :: ar_data
724 TYPE(arnoldi_control_type), POINTER :: control
725 TYPE(mp_comm_type) :: pcol_group
726
727 control => get_control(arnoldi_data)
728 pcol_group = control%pcol_group
729 local_comp = control%local_comp
730
731 ar_data => get_data_d(arnoldi_data)
732
733 ! create a local data copy to store the vectors and make Gram Schmidt a bit simpler
734 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
735 ALLOCATE (v_vec(nrow_local))
736 ALLOCATE (w_vec(nrow_local))
737 v_vec = 0.0_real_8; w_vec = 0.0_real_8
738 ar_data%Hessenberg = 0.0_real_8
739
740 IF (control%has_initial_vector) THEN
741 ! after calling the set routine the initial vector is stored in f_vec
742 CALL transfer_local_array_to_dbcsr_d (vectors%input_vec, ar_data%f_vec, nrow_local, control%local_comp)
743 ELSE
744 ! Setup the initial normalized random vector (sufficient if it only happens on proc_col 0)
745 CALL dbcsr_iterator_start(iter, vectors%input_vec)
746 DO WHILE (dbcsr_iterator_blocks_left(iter))
747 CALL dbcsr_iterator_next_block(iter, row, col, data_vec, row_size=row_size, col_size=col_size)
748 iseed(1) = 2; iseed(2) = mod(row, 4095); iseed(3) = mod(col, 4095); iseed(4) = 11
749 CALL dlarnv(2, iseed, row_size*col_size, data_vec)
750 END DO
751 CALL dbcsr_iterator_stop(iter)
752 END IF
753
754 CALL transfer_dbcsr_to_local_array_d (vectors%input_vec, v_vec, nrow_local, control%local_comp)
755
756 ! compute the vector norm of the reandom vectorm, get it real valued as well (rnorm)
757 CALL compute_norms_d (v_vec, norm, rnorm, control%pcol_group)
758
759 IF (rnorm == 0) rnorm = 1 ! catch case where this rank has no actual data
760 CALL dbcsr_scale(vectors%input_vec, real(1.0, real_8)/rnorm)
761
762 CALL dbcsr_matrix_colvec_multiply(matrix(1)%matrix, vectors%input_vec, vectors%result_vec, 1.0_real_8, &
763 0.0_real_8, vectors%rep_row_vec, vectors%rep_col_vec)
764
765 CALL transfer_dbcsr_to_local_array_d (vectors%result_vec, w_vec, nrow_local, control%local_comp)
766
767 ar_data%rho_scale = 0.0_real_8
768 ar_data%rho_scale = dot_product(v_vec, w_vec)
769 CALL pcol_group%sum(ar_data%rho_scale)
770
771 CALL dbcsr_matrix_colvec_multiply(matrix(2)%matrix, vectors%input_vec, vectors%result_vec, 1.0_real_8, &
772 0.0_real_8, vectors%rep_row_vec, vectors%rep_col_vec)
773
774 CALL transfer_dbcsr_to_local_array_d (vectors%result_vec, w_vec, nrow_local, control%local_comp)
775
776 denom = 0.0_real_8
777 denom = dot_product(v_vec, w_vec)
778 CALL pcol_group%sum(denom)
779 IF (control%myproc == 0) ar_data%rho_scale = ar_data%rho_scale/denom
780 CALL control%mp_group%bcast(ar_data%rho_scale, 0)
781
782 ! if the maximum ev is requested we need to optimize with -A-rho*B
783 CALL dbcsr_copy(matrix_arnoldi(1)%matrix, matrix(1)%matrix)
784 CALL dbcsr_add(matrix_arnoldi(1)%matrix, matrix(2)%matrix, 1.0_real_8, -ar_data%rho_scale)
785
786 ar_data%x_vec = v_vec
787
788 END SUBROUTINE gev_arnoldi_init_d
789
790! **************************************************************************************************
791!> \brief builds the basis rothogonal wrt. the metric.
792!> The structure looks similar to normal arnoldi but norms, vectors and
793!> matrix_vector products are very differently defined. Therefore it is
794!> cleaner to put it in a separate subroutine to avoid confusion
795!> \param matrix ...
796!> \param vectors ...
797!> \param arnoldi_data ...
798! **************************************************************************************************
799 SUBROUTINE gev_build_subspace_d (matrix, vectors, arnoldi_data)
800 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
801 TYPE(m_x_v_vectors_type) :: vectors
802 TYPE(arnoldi_data_type) :: arnoldi_data
803
804 INTEGER :: j, ncol_local, nrow_local
805 TYPE(arnoldi_control_type), POINTER :: control
806 TYPE(arnoldi_data_d_type), POINTER :: ar_data
807 REAL(kind=real_8) :: norm
808 REAL(kind=real_8), ALLOCATABLE, DIMENSION(:) :: h_vec, s_vec, v_vec, w_vec
809 REAL(kind=real_8), ALLOCATABLE, DIMENSION(:, :) :: zmat, czmat, bzmat
810 TYPE(mp_comm_type) :: pcol_group
811
812 ar_data => get_data_d(arnoldi_data)
813 control => get_control(arnoldi_data)
814 control%converged = .false.
815 pcol_group = control%pcol_group
816
817 ! create the vectors required during the iterations
818 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
819 ALLOCATE (v_vec(nrow_local)); ALLOCATE (w_vec(nrow_local))
820 v_vec = 0.0_real_8; w_vec = 0.0_real_8
821 ALLOCATE (s_vec(control%max_iter)); ALLOCATE (h_vec(control%max_iter))
822 ALLOCATE (zmat(nrow_local, control%max_iter)); ALLOCATE (czmat(nrow_local, control%max_iter))
823 ALLOCATE (bzmat(nrow_local, control%max_iter))
824
825 CALL transfer_local_array_to_dbcsr_d (vectors%input_vec, ar_data%x_vec, nrow_local, control%local_comp)
826 CALL dbcsr_matrix_colvec_multiply(matrix(2)%matrix, vectors%input_vec, vectors%result_vec, 1.0_real_8, &
827 0.0_real_8, vectors%rep_row_vec, vectors%rep_col_vec)
828 CALL transfer_dbcsr_to_local_array_d (vectors%result_vec, bzmat(:, 1), nrow_local, control%local_comp)
829
830 norm = 0.0_real_8
831 norm = dot_product(ar_data%x_vec, bzmat(:, 1))
832 CALL pcol_group%sum(norm)
833 IF (control%local_comp) THEN
834 zmat(:, 1) = ar_data%x_vec/sqrt(norm); bzmat(:, 1) = bzmat(:, 1)/sqrt(norm)
835 END IF
836
837 DO j = 1, control%max_iter
838 control%current_step = j
839 CALL transfer_local_array_to_dbcsr_d (vectors%input_vec, zmat(:, j), nrow_local, control%local_comp)
840 CALL dbcsr_matrix_colvec_multiply(matrix(1)%matrix, vectors%input_vec, vectors%result_vec, 1.0_real_8, &
841 0.0_real_8, vectors%rep_row_vec, vectors%rep_col_vec)
842 CALL transfer_dbcsr_to_local_array_d (vectors%result_vec, czmat(:, j), nrow_local, control%local_comp)
843 w_vec(:) = czmat(:, j)
844
845 ! Let's do the orthonormalization, to get the new f_vec. First try the Gram Schmidt scheme
846 CALL gram_schmidt_ortho_d (h_vec, ar_data%f_vec, s_vec, w_vec, nrow_local, j, &
847 bzmat, zmat, control%local_comp, control%pcol_group)
848
849 ! A bit more expensive but simpliy always top up with a DGKS correction, otherwise numerics
850 ! can becom a problem later on, there is probably a good check, but we don't perform it
851 CALL dgks_ortho_d (h_vec, ar_data%f_vec, s_vec, nrow_local, j, bzmat, &
852 zmat, control%local_comp, control%pcol_group)
853
854 CALL transfer_local_array_to_dbcsr_d (vectors%input_vec, ar_data%f_vec, nrow_local, control%local_comp)
855 CALL dbcsr_matrix_colvec_multiply(matrix(2)%matrix, vectors%input_vec, vectors%result_vec, 1.0_real_8, &
856 0.0_real_8, vectors%rep_row_vec, vectors%rep_col_vec)
857 CALL transfer_dbcsr_to_local_array_d (vectors%result_vec, v_vec, nrow_local, control%local_comp)
858 norm = 0.0_real_8
859 norm = dot_product(ar_data%f_vec, v_vec)
860 CALL pcol_group%sum(norm)
861
862 IF (control%myproc == 0) control%converged = real(norm, real_8) .LT. epsilon(real(1.0, real_8))
863 CALL control%mp_group%bcast(control%converged, 0)
864 IF (control%converged) EXIT
865 IF (j == control%max_iter - 1) EXIT
866
867 IF (control%local_comp) THEN
868 zmat(:, j + 1) = ar_data%f_vec/sqrt(norm); bzmat(:, j + 1) = v_vec(:)/sqrt(norm)
869 END IF
870 END DO
871
872! getting a bit more complicated here as the final matrix is again a product which has to be computed with the
873! distributed vectors, therefore a sum along the first proc_col is necessary. As we want that matrix everywhere,
874! we set it to zero before and compute the distributed product only on the first col and then sum over the full grid
875 ar_data%Hessenberg = 0.0_real_8
876 IF (control%local_comp) THEN
877 ar_data%Hessenberg(1:control%current_step, 1:control%current_step) = &
878 matmul(transpose(czmat(:, 1:control%current_step)), zmat(:, 1:control%current_step))
879 END IF
880 CALL control%mp_group%sum(ar_data%Hessenberg)
881
882 ar_data%local_history = zmat
883 ! broadcast the Hessenberg matrix so we don't need to care later on
884
885 DEALLOCATE (v_vec); DEALLOCATE (w_vec); DEALLOCATE (s_vec); DEALLOCATE (h_vec); DEALLOCATE (czmat);
886 DEALLOCATE (zmat); DEALLOCATE (bzmat)
887
888 END SUBROUTINE gev_build_subspace_d
889
890! **************************************************************************************************
891!> \brief Updates all data after an inner loop of the generalized ev arnoldi.
892!> Updates rho and C=A-rho*B accordingly.
893!> As an update scheme is used for he ev, the output ev has to be replaced
894!> with the updated one.
895!> Furthermore a convergence check is performed. The mv product could
896!> be skiiped by making clever use of the precomputed data, However,
897!> it is most likely not worth the effort.
898!> \param matrix ...
899!> \param matrix_arnoldi ...
900!> \param vectors ...
901!> \param arnoldi_data ...
902! **************************************************************************************************
903 SUBROUTINE gev_update_data_d (matrix, matrix_arnoldi, vectors, arnoldi_data)
904 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
905 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix_arnoldi
906 TYPE(m_x_v_vectors_type) :: vectors
907 TYPE(arnoldi_data_type) :: arnoldi_data
908
909 TYPE(arnoldi_control_type), POINTER :: control
910 TYPE(arnoldi_data_d_type), POINTER :: ar_data
911 INTEGER :: ncol_local, nrow_local, ind, i
912 REAL(kind=real_8), ALLOCATABLE, DIMENSION(:) :: v_vec
913 REAL(real_8) :: rnorm
914 REAL(kind=real_8) :: norm
915 COMPLEX(real_8) :: val
916
917 control => get_control(arnoldi_data)
918
919 ar_data => get_data_d(arnoldi_data)
920
921! compute the new shift, hack around the problem templating the conversion
922 val = ar_data%evals(control%selected_ind(1))
923 ar_data%rho_scale = ar_data%rho_scale + real(val,real_8)
924! compute the new eigenvector / initial guess for the next arnoldi loop
925 ar_data%x_vec = 0.0_real_8
926 DO i = 1, control%current_step
927 val = ar_data%revec(i, control%selected_ind(1))
928 ar_data%x_vec(:) = ar_data%x_vec(:) + ar_data%local_history(:, i)*real(val,real_8)
929 END DO
930! ar_data%x_vec(:)=MATMUL(ar_data%local_history(:,1:control%current_step),&
931! ar_data%revec(1:control%current_step,control%selected_ind(1)))
932
933! update the C-matrix (A-rho*B), if the maximum value is requested we have to use -A-rho*B
934 CALL dbcsr_copy(matrix_arnoldi(1)%matrix, matrix(1)%matrix)
935 CALL dbcsr_add(matrix_arnoldi(1)%matrix, matrix(2)%matrix, 1.0_real_8, -ar_data%rho_scale)
936
937! compute convergence and set the correct eigenvalue and eigenvector
938 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
939 IF (ncol_local > 0) THEN
940 ALLOCATE (v_vec(nrow_local))
941 CALL compute_norms_d (ar_data%x_vec, norm, rnorm, control%pcol_group)
942 v_vec(:) = ar_data%x_vec(:)/rnorm
943 END IF
944 CALL transfer_local_array_to_dbcsr_d (vectors%input_vec, v_vec, nrow_local, control%local_comp)
945 CALL dbcsr_matrix_colvec_multiply(matrix_arnoldi(1)%matrix, vectors%input_vec, vectors%result_vec, 1.0_real_8, &
946 0.0_real_8, vectors%rep_row_vec, vectors%rep_col_vec)
947 CALL transfer_dbcsr_to_local_array_d (vectors%result_vec, v_vec, nrow_local, control%local_comp)
948 IF (ncol_local > 0) THEN
949 CALL compute_norms_d (v_vec, norm, rnorm, control%pcol_group)
950 ! check convergence
951 control%converged = rnorm .LT. control%threshold
952 DEALLOCATE (v_vec)
953 END IF
954 ! and broadcast the real eigenvalue
955 CALL control%mp_group%bcast(control%converged, 0)
956 ind = control%selected_ind(1)
957 CALL control%mp_group%bcast(ar_data%rho_scale, 0)
958
959! Again the maximum value request is done on -A therefore the eigenvalue needs the opposite sign
960 ar_data%evals(ind) = ar_data%rho_scale
961
962 END SUBROUTINE gev_update_data_d
963! **************************************************************************************************
964!> \brief Call the correct eigensolver, in the arnoldi method only the right
965!> eigenvectors are used. Lefts are created here but dumped immediately
966!> \param arnoldi_data ...
967! **************************************************************************************************
968 SUBROUTINE compute_evals_z (arnoldi_data)
969 TYPE(arnoldi_data_type) :: arnoldi_data
970
971 COMPLEX(real_8), DIMENSION(:, :), ALLOCATABLE :: levec
972 TYPE(arnoldi_data_z_type), POINTER :: ar_data
973 INTEGER :: ndim
974 TYPE(arnoldi_control_type), POINTER :: control
975
976 ar_data => get_data_z(arnoldi_data)
977 control => get_control(arnoldi_data)
978 ndim = control%current_step
979 ALLOCATE (levec(ndim, ndim))
980
981! Needs antoher interface as the calls to real and complex geev differ (sucks!)
982! only perform the diagonalization on processors which hold data
983 IF (control%generalized_ev) THEN
984 CALL arnoldi_symm_local_diag('V', ar_data%Hessenberg(1:ndim, 1:ndim), ndim, &
985 ar_data%evals(1:ndim), ar_data%revec(1:ndim, 1:ndim))
986 ELSE
987 IF (control%symmetric) THEN
988 CALL arnoldi_tridiag_local_diag('N', 'V', ar_data%Hessenberg(1:ndim, 1:ndim), ndim, &
989 ar_data%evals(1:ndim), ar_data%revec(1:ndim, 1:ndim), levec)
990 ELSE
991 CALL arnoldi_general_local_diag('N', 'V', ar_data%Hessenberg(1:ndim, 1:ndim), ndim, &
992 ar_data%evals(1:ndim), ar_data%revec(1:ndim, 1:ndim), levec)
993 END IF
994 END IF
995
996 DEALLOCATE (levec)
997
998 END SUBROUTINE compute_evals_z
999
1000! **************************************************************************************************
1001!> \brief Initialization of the arnoldi vector. Here a random vector is used,
1002!> might be generalized in the future
1003!> \param matrix ...
1004!> \param vectors ...
1005!> \param arnoldi_data ...
1006! **************************************************************************************************
1007 SUBROUTINE arnoldi_init_z (matrix, vectors, arnoldi_data)
1008 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
1009 TYPE(m_x_v_vectors_type) :: vectors
1010 TYPE(arnoldi_data_type) :: arnoldi_data
1011
1012 INTEGER :: i, iseed(4), row_size, col_size, &
1013 nrow_local, ncol_local, &
1014 row, col
1015 REAL(real_8) :: rnorm
1016 TYPE(dbcsr_iterator_type) :: iter
1017 COMPLEX(kind=real_8) :: norm
1018 COMPLEX(kind=real_8), DIMENSION(:), ALLOCATABLE :: &
1019 v_vec, w_vec
1020 COMPLEX(kind=real_8), DIMENSION(:), POINTER :: data_vec
1021 LOGICAL :: local_comp
1022 TYPE(arnoldi_data_z_type), POINTER :: ar_data
1023 TYPE(arnoldi_control_type), POINTER :: control
1024 TYPE(mp_comm_type) :: pcol_group
1025
1026 control => get_control(arnoldi_data)
1027 pcol_group = control%pcol_group
1028 local_comp = control%local_comp
1029
1030 ar_data => get_data_z(arnoldi_data)
1031
1032 ! create a local data copy to store the vectors and make Gram Schmidt a bit simpler
1033 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
1034 ALLOCATE (v_vec(nrow_local))
1035 ALLOCATE (w_vec(nrow_local))
1036 v_vec = cmplx(0.0, 0.0, real_8); w_vec = cmplx(0.0, 0.0, real_8)
1037 ar_data%Hessenberg = cmplx(0.0, 0.0, real_8)
1038
1039 IF (control%has_initial_vector) THEN
1040 ! after calling the set routine the initial vector is stored in f_vec
1041 CALL transfer_local_array_to_dbcsr_z (vectors%input_vec, ar_data%f_vec, nrow_local, control%local_comp)
1042 ELSE
1043 ! Setup the initial normalized random vector (sufficient if it only happens on proc_col 0)
1044 CALL dbcsr_iterator_start(iter, vectors%input_vec)
1045 DO WHILE (dbcsr_iterator_blocks_left(iter))
1046 CALL dbcsr_iterator_next_block(iter, row, col, data_vec, row_size=row_size, col_size=col_size)
1047 iseed(1) = 2; iseed(2) = mod(row, 4095); iseed(3) = mod(col, 4095); iseed(4) = 11
1048 CALL zlarnv(2, iseed, row_size*col_size, data_vec)
1049 END DO
1050 CALL dbcsr_iterator_stop(iter)
1051 END IF
1052
1053 CALL transfer_dbcsr_to_local_array_z (vectors%input_vec, v_vec, nrow_local, control%local_comp)
1054
1055 ! compute the vector norm of the random vectorm, get it real valued as well (rnorm)
1056 CALL compute_norms_z (v_vec, norm, rnorm, control%pcol_group)
1057
1058 IF (rnorm == 0) rnorm = 1 ! catch case where this rank has no actual data
1059 CALL dbcsr_scale(vectors%input_vec, real(1.0, real_8)/rnorm)
1060
1061 ! Everything prepared, initialize the Arnoldi iteration
1062 CALL transfer_dbcsr_to_local_array_z (vectors%input_vec, v_vec, nrow_local, control%local_comp)
1063
1064 ! This permits to compute the subspace of a matrix which is a product of multiple matrices
1065 DO i = 1, SIZE(matrix)
1066 CALL dbcsr_matrix_colvec_multiply(matrix(i)%matrix, vectors%input_vec, vectors%result_vec, cmplx(1.0, 0.0, real_8), &
1067 cmplx(0.0, 0.0, real_8), vectors%rep_row_vec, vectors%rep_col_vec)
1068 CALL dbcsr_copy(vectors%input_vec, vectors%result_vec)
1069 END DO
1070
1071 CALL transfer_dbcsr_to_local_array_z (vectors%result_vec, w_vec, nrow_local, control%local_comp)
1072
1073 ! Put the projection into the Hessenberg matrix, and make the vectors orthonormal
1074 ar_data%Hessenberg(1, 1) = dot_product(v_vec, w_vec)
1075 CALL pcol_group%sum(ar_data%Hessenberg(1, 1))
1076 ar_data%f_vec = w_vec - v_vec*ar_data%Hessenberg(1, 1)
1077
1078 ar_data%local_history(:, 1) = v_vec(:)
1079
1080 ! We did the first step in here so we should set the current step for the subspace generation accordingly
1081 control%current_step = 1
1082
1083 DEALLOCATE (v_vec, w_vec)
1084
1085 END SUBROUTINE arnoldi_init_z
1086
1087! **************************************************************************************************
1088!> \brief Alogorithm for the implicit restarts in the arnoldi method
1089!> this is an early implementation which scales subspace size^4
1090!> by replacing the lapack calls with direct math the
1091!> QR and gemms can be made linear and a N^2 sacling will be acchieved
1092!> however this already sets the framework but should be used with care
1093!> \param arnoldi_data ...
1094! **************************************************************************************************
1095 SUBROUTINE arnoldi_iram_z (arnoldi_data)
1096 TYPE(arnoldi_data_type) :: arnoldi_data
1097
1098 TYPE(arnoldi_data_z_type), POINTER :: ar_data
1099 TYPE(arnoldi_control_type), POINTER :: control
1100 COMPLEX(real_8), DIMENSION(:, :), ALLOCATABLE :: tmp_mat, safe_mat, q, tmp_mat1
1101 COMPLEX(real_8), DIMENSION(:), ALLOCATABLE :: work, tau, work_measure
1102 INTEGER :: msize, lwork, i, j, info, nwant
1103 COMPLEX(kind=real_8) :: beta, sigma
1104 COMPLEX(kind=real_8), DIMENSION(:, :), ALLOCATABLE :: qdata
1105
1106 ! This is just a terribly inefficient implementation but I hope it is correct and might serve as a reference
1107 ar_data => get_data_z(arnoldi_data)
1108 control => get_control(arnoldi_data)
1109 msize = control%current_step
1110 nwant = control%nval_out
1111 ALLOCATE (tmp_mat(msize, msize)); ALLOCATE (safe_mat(msize, msize))
1112 ALLOCATE (q(msize, msize)); ALLOCATE (tmp_mat1(msize, msize))
1113 ALLOCATE (work_measure(1))
1114 ALLOCATE (tau(msize)); ALLOCATE (qdata(msize, msize))
1115 !make Q identity
1116 q = cmplx(0.0, 0.0, real_8)
1117 DO i = 1, msize
1118 q(i, i) = cmplx(1.0, 0.0, real_8)
1119 END DO
1120
1121 ! Looks a bit odd, but safe_mat will contain the result in the end, while tmpmat gets violated by lapack
1122 CALL convert_matrix(tmp_mat, ar_data%Hessenberg(1:msize, 1:msize))
1123 safe_mat(:, :) = tmp_mat(:, :)
1124
1125 DO i = 1, msize
1126 ! A bit a strange check but in the end we only want to shift the unwanted evals
1127 IF (any(control%selected_ind == i)) cycle
1128 ! Here we shift the matrix by subtracting unwanted evals from the diagonal
1129 DO j = 1, msize
1130 tmp_mat(j, j) = tmp_mat(j, j) - ar_data%evals(i)
1131 END DO
1132 ! Now we repair the damage by QR factorizing
1133 lwork = -1
1134 CALL zgeqrf(msize, msize, tmp_mat, msize, tau, work_measure, lwork, info)
1135 lwork = int(work_measure(1))
1136 IF (ALLOCATED(work)) THEN
1137 IF (SIZE(work) .LT. lwork) THEN
1138 DEALLOCATE (work)
1139 END IF
1140 END IF
1141 IF (.NOT. ALLOCATED(work)) ALLOCATE (work(lwork))
1142 CALL zgeqrf(msize, msize, tmp_mat, msize, tau, work, lwork, info)
1143 ! Ask Lapack to reconstruct Q from its own way of storing data (tmpmat will contain Q)
1144 CALL zungqr(msize, msize, msize, tmp_mat, msize, tau, work, lwork, info)
1145 ! update Q=Q*Q_current
1146 tmp_mat1(:, :) = q(:, :)
1147 CALL zgemm('N', 'N', msize, msize, msize, cmplx(1.0, 0.0, real_8), tmp_mat1, msize, tmp_mat, msize, cmplx(0.0, 0.0,&
1148 & real_8), &
1149 q, msize)
1150 ! Update H=(Q*)HQ
1151 CALL zgemm('C', 'N', msize, msize, msize, cmplx(1.0, 0.0, real_8), tmp_mat, msize, safe_mat, msize, cmplx(0.0, 0.0,&
1152 & real_8), &
1153 tmp_mat1, msize)
1154 CALL zgemm('N', 'N', msize, msize, msize, cmplx(1.0, 0.0, real_8), tmp_mat1, msize, tmp_mat, msize, cmplx(0.0, 0.0,&
1155 & real_8), &
1156 safe_mat, msize)
1157
1158 ! this one is crucial for numerics not to accumulate noise in the subdiagonals
1159 DO j = 1, msize
1160 safe_mat(j + 2:msize, j) = cmplx(0.0, 0.0, real_8)
1161 END DO
1162 tmp_mat(:, :) = safe_mat(:, :)
1163 END DO
1164
1165 ! Now we can compute our restart quantities
1166 ar_data%Hessenberg = cmplx(0.0, 0.0, real_8)
1167 CALL convert_matrix(ar_data%Hessenberg(1:msize, 1:msize), safe_mat)
1168 CALL convert_matrix(qdata, q)
1169
1170 beta = ar_data%Hessenberg(nwant + 1, nwant); sigma = qdata(msize, nwant)
1171
1172 !update the residuum and the basis vectors
1173 IF (control%local_comp) THEN
1174 ar_data%f_vec = matmul(ar_data%local_history(:, 1:msize), qdata(1:msize, nwant + 1))*beta + ar_data%f_vec(:)*sigma
1175 ar_data%local_history(:, 1:nwant) = matmul(ar_data%local_history(:, 1:msize), qdata(1:msize, 1:nwant))
1176 END IF
1177 ! Set the current step to nwant so the subspace build knows where to start
1178 control%current_step = nwant
1179
1180 DEALLOCATE (tmp_mat, safe_mat, q, qdata, tmp_mat1, work, tau, work_measure)
1181
1182 END SUBROUTINE arnoldi_iram_z
1183
1184! **************************************************************************************************
1185!> \brief Here we create the Krylov subspace and fill the Hessenberg matrix
1186!> convergence check is only performed on subspace convergence
1187!> Gram Schidt is used to orthonogonalize.
1188!> If this is numericall not sufficient a Daniel, Gragg, Kaufman and Steward
1189!> correction is performed
1190!> \param matrix ...
1191!> \param vectors ...
1192!> \param arnoldi_data ...
1193! **************************************************************************************************
1194 SUBROUTINE build_subspace_z (matrix, vectors, arnoldi_data)
1195 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
1196 TYPE(m_x_v_vectors_type), TARGET :: vectors
1197 TYPE(arnoldi_data_type) :: arnoldi_data
1198
1199 INTEGER :: i, j, ncol_local, nrow_local
1200 REAL(real_8) :: rnorm
1201 TYPE(arnoldi_control_type), POINTER :: control
1202 TYPE(arnoldi_data_z_type), POINTER :: ar_data
1203 COMPLEX(kind=real_8) :: norm
1204 COMPLEX(kind=real_8), ALLOCATABLE, DIMENSION(:) :: h_vec, s_vec, v_vec, w_vec
1205 TYPE(dbcsr_type), POINTER :: input_vec, result_vec, swap_vec
1206
1207 ar_data => get_data_z(arnoldi_data)
1208 control => get_control(arnoldi_data)
1209 control%converged = .false.
1210
1211 ! create the vectors required during the iterations
1212 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
1213 ALLOCATE (v_vec(nrow_local)); ALLOCATE (w_vec(nrow_local))
1214 v_vec = cmplx(0.0, 0.0, real_8); w_vec = cmplx(0.0, 0.0, real_8)
1215 ALLOCATE (s_vec(control%max_iter)); ALLOCATE (h_vec(control%max_iter))
1216
1217 DO j = control%current_step, control%max_iter - 1
1218
1219 ! compute the vector norm of the residuum, get it real valued as well (rnorm)
1220 CALL compute_norms_z (ar_data%f_vec, norm, rnorm, control%pcol_group)
1221
1222 ! check convergence and inform everybody about it, a bit annoying to talk to everybody because of that
1223 IF (control%myproc == 0) control%converged = rnorm .LT. real(control%threshold, real_8)
1224 CALL control%mp_group%bcast(control%converged, 0)
1225 IF (control%converged) EXIT
1226
1227 ! transfer normalized residdum to history and its norm to the Hessenberg matrix
1228 IF (rnorm == 0) rnorm = 1 ! catch case where this rank has no actual data
1229 v_vec(:) = ar_data%f_vec(:)/rnorm; ar_data%local_history(:, j + 1) = v_vec(:); ar_data%Hessenberg(j + 1, j) = norm
1230
1231 input_vec => vectors%input_vec
1232 result_vec => vectors%result_vec
1233 CALL transfer_local_array_to_dbcsr_z (input_vec, v_vec, nrow_local, control%local_comp)
1234
1235 ! This permits to compute the subspace of a matrix which is a product of two matrices
1236 DO i = 1, SIZE(matrix)
1237 CALL dbcsr_matrix_colvec_multiply(matrix(i)%matrix, input_vec, result_vec, cmplx(1.0, 0.0, real_8), &
1238 cmplx(0.0, 0.0, real_8), vectors%rep_row_vec, vectors%rep_col_vec)
1239 swap_vec => input_vec
1240 input_vec => result_vec
1241 result_vec => swap_vec
1242 END DO
1243
1244 CALL transfer_dbcsr_to_local_array_z (input_vec, w_vec, nrow_local, control%local_comp)
1245
1246 ! Let's do the orthonormalization, to get the new f_vec. First try the Gram Schmidt scheme
1247 CALL gram_schmidt_ortho_z (h_vec, ar_data%f_vec, s_vec, w_vec, nrow_local, j + 1, &
1248 ar_data%local_history, ar_data%local_history, control%local_comp, control%pcol_group)
1249
1250 ! A bit more expensive but simply always top up with a DGKS correction, otherwise numerics
1251 ! can become a problem later on, there is probably a good check whether it's necessary, but we don't perform it
1252 CALL dgks_ortho_z (h_vec, ar_data%f_vec, s_vec, nrow_local, j + 1, ar_data%local_history, &
1253 ar_data%local_history, control%local_comp, control%pcol_group)
1254 ! Finally we can put the projections into our Hessenberg matrix
1255 ar_data%Hessenberg(1:j + 1, j + 1) = h_vec(1:j + 1)
1256 control%current_step = j + 1
1257 END DO
1258
1259 ! compute the vector norm of the final residuum and put it in to Hessenberg
1260 CALL compute_norms_z (ar_data%f_vec, norm, rnorm, control%pcol_group)
1261 ar_data%Hessenberg(control%current_step + 1, control%current_step) = norm
1262
1263 ! broadcast the Hessenberg matrix so we don't need to care later on
1264 CALL control%mp_group%bcast(ar_data%Hessenberg, 0)
1265
1266 DEALLOCATE (v_vec, w_vec, h_vec, s_vec)
1267
1268 END SUBROUTINE build_subspace_z
1269
1270! **************************************************************************************************
1271!> \brief Helper routine to transfer the all data of a dbcsr matrix to a local array
1272!> \param vec ...
1273!> \param array ...
1274!> \param n ...
1275!> \param is_local ...
1276! **************************************************************************************************
1277 SUBROUTINE transfer_dbcsr_to_local_array_z (vec, array, n, is_local)
1278 TYPE(dbcsr_type) :: vec
1279 COMPLEX(kind=real_8), DIMENSION(:) :: array
1280 INTEGER :: n
1281 LOGICAL :: is_local
1282 COMPLEX(kind=real_8), DIMENSION(:), POINTER :: data_vec
1283
1284 data_vec => dbcsr_get_data_p(vec, select_data_type=cmplx(0.0, 0.0, real_8))
1285 IF (is_local) array(1:n) = data_vec(1:n)
1286
1287 END SUBROUTINE transfer_dbcsr_to_local_array_z
1288
1289! **************************************************************************************************
1290!> \brief The inverse routine transferring data back from an array to a dbcsr
1291!> \param vec ...
1292!> \param array ...
1293!> \param n ...
1294!> \param is_local ...
1295! **************************************************************************************************
1296 SUBROUTINE transfer_local_array_to_dbcsr_z (vec, array, n, is_local)
1297 TYPE(dbcsr_type) :: vec
1298 COMPLEX(kind=real_8), DIMENSION(:) :: array
1299 INTEGER :: n
1300 LOGICAL :: is_local
1301 COMPLEX(kind=real_8), DIMENSION(:), POINTER :: data_vec
1302
1303 data_vec => dbcsr_get_data_p(vec, select_data_type=cmplx(0.0, 0.0, real_8))
1304 IF (is_local) data_vec(1:n) = array(1:n)
1305
1306 END SUBROUTINE transfer_local_array_to_dbcsr_z
1307
1308! **************************************************************************************************
1309!> \brief Gram-Schmidt in matrix vector form
1310!> \param h_vec ...
1311!> \param f_vec ...
1312!> \param s_vec ...
1313!> \param w_vec ...
1314!> \param nrow_local ...
1315!> \param j ...
1316!> \param local_history ...
1317!> \param reorth_mat ...
1318!> \param local_comp ...
1319!> \param pcol_group ...
1320! **************************************************************************************************
1321 SUBROUTINE gram_schmidt_ortho_z (h_vec, f_vec, s_vec, w_vec, nrow_local, &
1322 j, local_history, reorth_mat, local_comp, pcol_group)
1323 COMPLEX(kind=real_8), DIMENSION(:) :: h_vec, s_vec, f_vec, w_vec
1324 COMPLEX(kind=real_8), DIMENSION(:, :) :: local_history, reorth_mat
1325 INTEGER :: nrow_local, j
1326 TYPE(mp_comm_type), INTENT(IN) :: pcol_group
1327 LOGICAL :: local_comp
1328
1329 CHARACTER(LEN=*), PARAMETER :: routinen = 'Gram_Schmidt_ortho_z'
1330 INTEGER :: handle
1331
1332 CALL timeset(routinen, handle)
1333
1334 ! Let's do the orthonormalization, first try the Gram Schmidt scheme
1335 h_vec = cmplx(0.0, 0.0, real_8); f_vec = cmplx(0.0, 0.0, real_8); s_vec = cmplx(0.0, 0.0, real_8)
1336 IF (local_comp) CALL zgemv('T', nrow_local, j, cmplx(1.0, 0.0, real_8), local_history, &
1337 nrow_local, w_vec, 1, cmplx(0.0, 0.0, real_8), h_vec, 1)
1338 CALL pcol_group%sum(h_vec(1:j))
1339
1340 IF (local_comp) CALL zgemv('N', nrow_local, j, cmplx(1.0, 0.0, real_8), reorth_mat, &
1341 nrow_local, h_vec, 1, cmplx(0.0, 0.0, real_8), f_vec, 1)
1342 f_vec(:) = w_vec(:) - f_vec(:)
1343
1344 CALL timestop(handle)
1345
1346 END SUBROUTINE gram_schmidt_ortho_z
1347
1348! **************************************************************************************************
1349!> \brief Compute the Daniel, Gragg, Kaufman and Steward correction
1350!> \param h_vec ...
1351!> \param f_vec ...
1352!> \param s_vec ...
1353!> \param nrow_local ...
1354!> \param j ...
1355!> \param local_history ...
1356!> \param reorth_mat ...
1357!> \param local_comp ...
1358!> \param pcol_group ...
1359! **************************************************************************************************
1360 SUBROUTINE dgks_ortho_z (h_vec, f_vec, s_vec, nrow_local, j, &
1361 local_history, reorth_mat, local_comp, pcol_group)
1362 COMPLEX(kind=real_8), DIMENSION(:) :: h_vec, s_vec, f_vec
1363 COMPLEX(kind=real_8), DIMENSION(:, :) :: local_history, reorth_mat
1364 INTEGER :: nrow_local, j
1365 TYPE(mp_comm_type), INTENT(IN) :: pcol_group
1366
1367 CHARACTER(LEN=*), PARAMETER :: routinen = 'DGKS_ortho_z'
1368 LOGICAL :: local_comp
1369 INTEGER :: handle
1370
1371 CALL timeset(routinen, handle)
1372
1373 IF (local_comp) CALL zgemv('T', nrow_local, j, cmplx(1.0, 0.0, real_8), local_history, &
1374 nrow_local, f_vec, 1, cmplx(0.0, 0.0, real_8), s_vec, 1)
1375 CALL pcol_group%sum(s_vec(1:j))
1376 IF (local_comp) CALL zgemv('N', nrow_local, j, cmplx(-1.0, 0.0, real_8), reorth_mat, &
1377 nrow_local, s_vec, 1, cmplx(1.0, 0.0, real_8), f_vec, 1)
1378 h_vec(1:j) = h_vec(1:j) + s_vec(1:j)
1379
1380 CALL timestop(handle)
1381
1382 END SUBROUTINE dgks_ortho_z
1383
1384! **************************************************************************************************
1385!> \brief Compute the norm of a vector distributed along proc_col
1386!> as local arrays. Always return the real part next to the complex rep.
1387!> \param vec ...
1388!> \param norm ...
1389!> \param rnorm ...
1390!> \param pcol_group ...
1391! **************************************************************************************************
1392 SUBROUTINE compute_norms_z (vec, norm, rnorm, pcol_group)
1393 COMPLEX(kind=real_8), DIMENSION(:) :: vec
1394 REAL(real_8) :: rnorm
1395 COMPLEX(kind=real_8) :: norm
1396 TYPE(mp_comm_type), INTENT(IN) :: pcol_group
1397
1398 ! the norm is computed along the processor column
1399 norm = dot_product(vec, vec)
1400 CALL pcol_group%sum(norm)
1401 rnorm = sqrt(real(norm, real_8))
1402 norm=cmplx(rnorm, 0.0, real_8)
1403
1404 END SUBROUTINE compute_norms_z
1405
1406! **************************************************************************************************
1407!> \brief Computes the initial guess for the solution of the generalized eigenvalue
1408!> using the arnoldi method
1409!> \param matrix ...
1410!> \param matrix_arnoldi ...
1411!> \param vectors ...
1412!> \param arnoldi_data ...
1413! **************************************************************************************************
1414
1415 SUBROUTINE gev_arnoldi_init_z (matrix, matrix_arnoldi, vectors, arnoldi_data)
1416 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
1417 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix_arnoldi
1418 TYPE(m_x_v_vectors_type) :: vectors
1419 TYPE(arnoldi_data_type) :: arnoldi_data
1420
1421 INTEGER :: iseed(4), row_size, col_size, &
1422 nrow_local, ncol_local, &
1423 row, col
1424 REAL(real_8) :: rnorm
1425 TYPE(dbcsr_iterator_type) :: iter
1426 COMPLEX(kind=real_8) :: norm, denom
1427 COMPLEX(kind=real_8), DIMENSION(:), ALLOCATABLE :: &
1428 v_vec, w_vec
1429 COMPLEX(kind=real_8), DIMENSION(:), POINTER :: data_vec
1430 LOGICAL :: local_comp
1431 TYPE(arnoldi_data_z_type), POINTER :: ar_data
1432 TYPE(arnoldi_control_type), POINTER :: control
1433 TYPE(mp_comm_type) :: pcol_group
1434
1435 control => get_control(arnoldi_data)
1436 pcol_group = control%pcol_group
1437 local_comp = control%local_comp
1438
1439 ar_data => get_data_z(arnoldi_data)
1440
1441 ! create a local data copy to store the vectors and make Gram Schmidt a bit simpler
1442 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
1443 ALLOCATE (v_vec(nrow_local))
1444 ALLOCATE (w_vec(nrow_local))
1445 v_vec = cmplx(0.0, 0.0, real_8); w_vec = cmplx(0.0, 0.0, real_8)
1446 ar_data%Hessenberg = cmplx(0.0, 0.0, real_8)
1447
1448 IF (control%has_initial_vector) THEN
1449 ! after calling the set routine the initial vector is stored in f_vec
1450 CALL transfer_local_array_to_dbcsr_z (vectors%input_vec, ar_data%f_vec, nrow_local, control%local_comp)
1451 ELSE
1452 ! Setup the initial normalized random vector (sufficient if it only happens on proc_col 0)
1453 CALL dbcsr_iterator_start(iter, vectors%input_vec)
1454 DO WHILE (dbcsr_iterator_blocks_left(iter))
1455 CALL dbcsr_iterator_next_block(iter, row, col, data_vec, row_size=row_size, col_size=col_size)
1456 iseed(1) = 2; iseed(2) = mod(row, 4095); iseed(3) = mod(col, 4095); iseed(4) = 11
1457 CALL zlarnv(2, iseed, row_size*col_size, data_vec)
1458 END DO
1459 CALL dbcsr_iterator_stop(iter)
1460 END IF
1461
1462 CALL transfer_dbcsr_to_local_array_z (vectors%input_vec, v_vec, nrow_local, control%local_comp)
1463
1464 ! compute the vector norm of the reandom vectorm, get it real valued as well (rnorm)
1465 CALL compute_norms_z (v_vec, norm, rnorm, control%pcol_group)
1466
1467 IF (rnorm == 0) rnorm = 1 ! catch case where this rank has no actual data
1468 CALL dbcsr_scale(vectors%input_vec, real(1.0, real_8)/rnorm)
1469
1470 CALL dbcsr_matrix_colvec_multiply(matrix(1)%matrix, vectors%input_vec, vectors%result_vec, cmplx(1.0, 0.0, real_8), &
1471 cmplx(0.0, 0.0, real_8), vectors%rep_row_vec, vectors%rep_col_vec)
1472
1473 CALL transfer_dbcsr_to_local_array_z (vectors%result_vec, w_vec, nrow_local, control%local_comp)
1474
1475 ar_data%rho_scale = cmplx(0.0, 0.0, real_8)
1476 ar_data%rho_scale = dot_product(v_vec, w_vec)
1477 CALL pcol_group%sum(ar_data%rho_scale)
1478
1479 CALL dbcsr_matrix_colvec_multiply(matrix(2)%matrix, vectors%input_vec, vectors%result_vec, cmplx(1.0, 0.0, real_8), &
1480 cmplx(0.0, 0.0, real_8), vectors%rep_row_vec, vectors%rep_col_vec)
1481
1482 CALL transfer_dbcsr_to_local_array_z (vectors%result_vec, w_vec, nrow_local, control%local_comp)
1483
1484 denom = cmplx(0.0, 0.0, real_8)
1485 denom = dot_product(v_vec, w_vec)
1486 CALL pcol_group%sum(denom)
1487 IF (control%myproc == 0) ar_data%rho_scale = ar_data%rho_scale/denom
1488 CALL control%mp_group%bcast(ar_data%rho_scale, 0)
1489
1490 ! if the maximum ev is requested we need to optimize with -A-rho*B
1491 CALL dbcsr_copy(matrix_arnoldi(1)%matrix, matrix(1)%matrix)
1492 CALL dbcsr_add(matrix_arnoldi(1)%matrix, matrix(2)%matrix, cmplx(1.0, 0.0, real_8), -ar_data%rho_scale)
1493
1494 ar_data%x_vec = v_vec
1495
1496 END SUBROUTINE gev_arnoldi_init_z
1497
1498! **************************************************************************************************
1499!> \brief builds the basis rothogonal wrt. the metric.
1500!> The structure looks similar to normal arnoldi but norms, vectors and
1501!> matrix_vector products are very differently defined. Therefore it is
1502!> cleaner to put it in a separate subroutine to avoid confusion
1503!> \param matrix ...
1504!> \param vectors ...
1505!> \param arnoldi_data ...
1506! **************************************************************************************************
1507 SUBROUTINE gev_build_subspace_z (matrix, vectors, arnoldi_data)
1508 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
1509 TYPE(m_x_v_vectors_type) :: vectors
1510 TYPE(arnoldi_data_type) :: arnoldi_data
1511
1512 INTEGER :: j, ncol_local, nrow_local
1513 TYPE(arnoldi_control_type), POINTER :: control
1514 TYPE(arnoldi_data_z_type), POINTER :: ar_data
1515 COMPLEX(kind=real_8) :: norm
1516 COMPLEX(kind=real_8), ALLOCATABLE, DIMENSION(:) :: h_vec, s_vec, v_vec, w_vec
1517 COMPLEX(kind=real_8), ALLOCATABLE, DIMENSION(:, :) :: zmat, czmat, bzmat
1518 TYPE(mp_comm_type) :: pcol_group
1519
1520 ar_data => get_data_z(arnoldi_data)
1521 control => get_control(arnoldi_data)
1522 control%converged = .false.
1523 pcol_group = control%pcol_group
1524
1525 ! create the vectors required during the iterations
1526 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
1527 ALLOCATE (v_vec(nrow_local)); ALLOCATE (w_vec(nrow_local))
1528 v_vec = cmplx(0.0, 0.0, real_8); w_vec = cmplx(0.0, 0.0, real_8)
1529 ALLOCATE (s_vec(control%max_iter)); ALLOCATE (h_vec(control%max_iter))
1530 ALLOCATE (zmat(nrow_local, control%max_iter)); ALLOCATE (czmat(nrow_local, control%max_iter))
1531 ALLOCATE (bzmat(nrow_local, control%max_iter))
1532
1533 CALL transfer_local_array_to_dbcsr_z (vectors%input_vec, ar_data%x_vec, nrow_local, control%local_comp)
1534 CALL dbcsr_matrix_colvec_multiply(matrix(2)%matrix, vectors%input_vec, vectors%result_vec, cmplx(1.0, 0.0, real_8), &
1535 cmplx(0.0, 0.0, real_8), vectors%rep_row_vec, vectors%rep_col_vec)
1536 CALL transfer_dbcsr_to_local_array_z (vectors%result_vec, bzmat(:, 1), nrow_local, control%local_comp)
1537
1538 norm = cmplx(0.0, 0.0, real_8)
1539 norm = dot_product(ar_data%x_vec, bzmat(:, 1))
1540 CALL pcol_group%sum(norm)
1541 IF (control%local_comp) THEN
1542 zmat(:, 1) = ar_data%x_vec/sqrt(norm); bzmat(:, 1) = bzmat(:, 1)/sqrt(norm)
1543 END IF
1544
1545 DO j = 1, control%max_iter
1546 control%current_step = j
1547 CALL transfer_local_array_to_dbcsr_z (vectors%input_vec, zmat(:, j), nrow_local, control%local_comp)
1548 CALL dbcsr_matrix_colvec_multiply(matrix(1)%matrix, vectors%input_vec, vectors%result_vec, cmplx(1.0, 0.0, real_8), &
1549 cmplx(0.0, 0.0, real_8), vectors%rep_row_vec, vectors%rep_col_vec)
1550 CALL transfer_dbcsr_to_local_array_z (vectors%result_vec, czmat(:, j), nrow_local, control%local_comp)
1551 w_vec(:) = czmat(:, j)
1552
1553 ! Let's do the orthonormalization, to get the new f_vec. First try the Gram Schmidt scheme
1554 CALL gram_schmidt_ortho_z (h_vec, ar_data%f_vec, s_vec, w_vec, nrow_local, j, &
1555 bzmat, zmat, control%local_comp, control%pcol_group)
1556
1557 ! A bit more expensive but simpliy always top up with a DGKS correction, otherwise numerics
1558 ! can becom a problem later on, there is probably a good check, but we don't perform it
1559 CALL dgks_ortho_z (h_vec, ar_data%f_vec, s_vec, nrow_local, j, bzmat, &
1560 zmat, control%local_comp, control%pcol_group)
1561
1562 CALL transfer_local_array_to_dbcsr_z (vectors%input_vec, ar_data%f_vec, nrow_local, control%local_comp)
1563 CALL dbcsr_matrix_colvec_multiply(matrix(2)%matrix, vectors%input_vec, vectors%result_vec, cmplx(1.0, 0.0, real_8), &
1564 cmplx(0.0, 0.0, real_8), vectors%rep_row_vec, vectors%rep_col_vec)
1565 CALL transfer_dbcsr_to_local_array_z (vectors%result_vec, v_vec, nrow_local, control%local_comp)
1566 norm = cmplx(0.0, 0.0, real_8)
1567 norm = dot_product(ar_data%f_vec, v_vec)
1568 CALL pcol_group%sum(norm)
1569
1570 IF (control%myproc == 0) control%converged = real(norm, real_8) .LT. epsilon(real(1.0, real_8))
1571 CALL control%mp_group%bcast(control%converged, 0)
1572 IF (control%converged) EXIT
1573 IF (j == control%max_iter - 1) EXIT
1574
1575 IF (control%local_comp) THEN
1576 zmat(:, j + 1) = ar_data%f_vec/sqrt(norm); bzmat(:, j + 1) = v_vec(:)/sqrt(norm)
1577 END IF
1578 END DO
1579
1580! getting a bit more complicated here as the final matrix is again a product which has to be computed with the
1581! distributed vectors, therefore a sum along the first proc_col is necessary. As we want that matrix everywhere,
1582! we set it to zero before and compute the distributed product only on the first col and then sum over the full grid
1583 ar_data%Hessenberg = cmplx(0.0, 0.0, real_8)
1584 IF (control%local_comp) THEN
1585 ar_data%Hessenberg(1:control%current_step, 1:control%current_step) = &
1586 matmul(transpose(czmat(:, 1:control%current_step)), zmat(:, 1:control%current_step))
1587 END IF
1588 CALL control%mp_group%sum(ar_data%Hessenberg)
1589
1590 ar_data%local_history = zmat
1591 ! broadcast the Hessenberg matrix so we don't need to care later on
1592
1593 DEALLOCATE (v_vec); DEALLOCATE (w_vec); DEALLOCATE (s_vec); DEALLOCATE (h_vec); DEALLOCATE (czmat);
1594 DEALLOCATE (zmat); DEALLOCATE (bzmat)
1595
1596 END SUBROUTINE gev_build_subspace_z
1597
1598! **************************************************************************************************
1599!> \brief Updates all data after an inner loop of the generalized ev arnoldi.
1600!> Updates rho and C=A-rho*B accordingly.
1601!> As an update scheme is used for he ev, the output ev has to be replaced
1602!> with the updated one.
1603!> Furthermore a convergence check is performed. The mv product could
1604!> be skiiped by making clever use of the precomputed data, However,
1605!> it is most likely not worth the effort.
1606!> \param matrix ...
1607!> \param matrix_arnoldi ...
1608!> \param vectors ...
1609!> \param arnoldi_data ...
1610! **************************************************************************************************
1611 SUBROUTINE gev_update_data_z (matrix, matrix_arnoldi, vectors, arnoldi_data)
1612 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
1613 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix_arnoldi
1614 TYPE(m_x_v_vectors_type) :: vectors
1615 TYPE(arnoldi_data_type) :: arnoldi_data
1616
1617 TYPE(arnoldi_control_type), POINTER :: control
1618 TYPE(arnoldi_data_z_type), POINTER :: ar_data
1619 INTEGER :: ncol_local, nrow_local, ind, i
1620 COMPLEX(kind=real_8), ALLOCATABLE, DIMENSION(:) :: v_vec
1621 REAL(real_8) :: rnorm
1622 COMPLEX(kind=real_8) :: norm
1623 COMPLEX(real_8) :: val
1624
1625 control => get_control(arnoldi_data)
1626
1627 ar_data => get_data_z(arnoldi_data)
1628
1629! compute the new shift, hack around the problem templating the conversion
1630 val = ar_data%evals(control%selected_ind(1))
1631 ar_data%rho_scale = ar_data%rho_scale + val
1632! compute the new eigenvector / initial guess for the next arnoldi loop
1633 ar_data%x_vec = cmplx(0.0, 0.0, real_8)
1634 DO i = 1, control%current_step
1635 val = ar_data%revec(i, control%selected_ind(1))
1636 ar_data%x_vec(:) = ar_data%x_vec(:) + ar_data%local_history(:, i)*val
1637 END DO
1638! ar_data%x_vec(:)=MATMUL(ar_data%local_history(:,1:control%current_step),&
1639! ar_data%revec(1:control%current_step,control%selected_ind(1)))
1640
1641! update the C-matrix (A-rho*B), if the maximum value is requested we have to use -A-rho*B
1642 CALL dbcsr_copy(matrix_arnoldi(1)%matrix, matrix(1)%matrix)
1643 CALL dbcsr_add(matrix_arnoldi(1)%matrix, matrix(2)%matrix, cmplx(1.0, 0.0, real_8), -ar_data%rho_scale)
1644
1645! compute convergence and set the correct eigenvalue and eigenvector
1646 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
1647 IF (ncol_local > 0) THEN
1648 ALLOCATE (v_vec(nrow_local))
1649 CALL compute_norms_z (ar_data%x_vec, norm, rnorm, control%pcol_group)
1650 v_vec(:) = ar_data%x_vec(:)/rnorm
1651 END IF
1652 CALL transfer_local_array_to_dbcsr_z (vectors%input_vec, v_vec, nrow_local, control%local_comp)
1653 CALL dbcsr_matrix_colvec_multiply(matrix_arnoldi(1)%matrix, vectors%input_vec, vectors%result_vec, cmplx(1.0, 0.0,&
1654 & real_8), &
1655 cmplx(0.0, 0.0, real_8), vectors%rep_row_vec, vectors%rep_col_vec)
1656 CALL transfer_dbcsr_to_local_array_z (vectors%result_vec, v_vec, nrow_local, control%local_comp)
1657 IF (ncol_local > 0) THEN
1658 CALL compute_norms_z (v_vec, norm, rnorm, control%pcol_group)
1659 ! check convergence
1660 control%converged = rnorm .LT. control%threshold
1661 DEALLOCATE (v_vec)
1662 END IF
1663 ! and broadcast the real eigenvalue
1664 CALL control%mp_group%bcast(control%converged, 0)
1665 ind = control%selected_ind(1)
1666 CALL control%mp_group%bcast(ar_data%rho_scale, 0)
1667
1668! Again the maximum value request is done on -A therefore the eigenvalue needs the opposite sign
1669 ar_data%evals(ind) = ar_data%rho_scale
1670
1671 END SUBROUTINE gev_update_data_z
1672
1673END MODULE arnoldi_methods
arnoldi_geev::arnoldi_general_local_diag
Definition arnoldi_geev.F:28
arnoldi_geev::arnoldi_symm_local_diag
Definition arnoldi_geev.F:38
arnoldi_geev::arnoldi_tridiag_local_diag
Definition arnoldi_geev.F:33
dbcsr_vector::dbcsr_matrix_colvec_multiply
Definition dbcsr_vector.F:92
arnoldi_geev
provides a unified interface to lapack geev routines
Definition arnoldi_geev.F:15
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_data_c
type(arnoldi_data_c_type) function, pointer, public get_data_c(ar_data)
...
Definition arnoldi_types.F:180
arnoldi_types::get_data_s
type(arnoldi_data_s_type) function, pointer, public get_data_s(ar_data)
...
Definition arnoldi_types.F:154
arnoldi_types::get_data_z
type(arnoldi_data_z_type) function, pointer, public get_data_z(ar_data)
...
Definition arnoldi_types.F:167
arnoldi_types::has_d_real
logical function, public has_d_real(ar_data)
...
Definition arnoldi_types.F:258
arnoldi_types::get_data_d
type(arnoldi_data_d_type) function, pointer, public get_data_d(ar_data)
...
Definition arnoldi_types.F:141
arnoldi_types::get_control
type(arnoldi_control_type) function, pointer, public get_control(ar_data)
...
Definition arnoldi_types.F:245
arnoldi_types::has_d_cmplx
elemental logical function, public has_d_cmplx(ar_data)
...
Definition arnoldi_types.F:284
dbcsr_vector
operations for skinny matrices/vectors expressed in dbcsr form
Definition dbcsr_vector.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::real_8
integer, parameter, public real_8
Definition kinds.F:41
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_c_type
Definition arnoldi_types.F:74
arnoldi_types::arnoldi_data_d_type
Definition arnoldi_types.F:44
arnoldi_types::arnoldi_data_s_type
Definition arnoldi_types.F:54
arnoldi_types::arnoldi_data_type
Definition arnoldi_types.F:84
arnoldi_types::arnoldi_data_z_type
Definition arnoldi_types.F:64
arnoldi_types::m_x_v_vectors_type
Definition arnoldi_types.F:92
message_passing::mp_comm_type
Definition message_passing.F:189