(git:98357aa)
Loading...
Searching...
No Matches
iterate_matrix.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7! **************************************************************************************************
8!> \brief Routines useful for iterative matrix calculations
9!> \par History
10!> 2010.10 created [Joost VandeVondele]
11!> \author Joost VandeVondele
12! **************************************************************************************************
13MODULE iterate_matrix
14 USE arnoldi_api, ONLY: arnoldi_env_type,&
15 arnoldi_extremal
16 USE bibliography, ONLY: richters2018,&
17 cite_reference
18 USE cp_dbcsr_api, ONLY: &
19 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_get, &
20 dbcsr_distribution_type, dbcsr_filter, dbcsr_get_info, dbcsr_get_matrix_type, &
21 dbcsr_get_occupation, dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, &
22 dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry
23 USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
24 dbcsr_frobenius_norm,&
25 dbcsr_gershgorin_norm,&
26 dbcsr_get_diag,&
27 dbcsr_maxabs,&
28 dbcsr_set_diag,&
29 dbcsr_trace
30 USE cp_log_handling, ONLY: cp_get_default_logger,&
31 cp_logger_get_default_unit_nr,&
32 cp_logger_type
33 USE input_constants, ONLY: ls_scf_submatrix_sign_direct,&
34 ls_scf_submatrix_sign_direct_muadj,&
35 ls_scf_submatrix_sign_direct_muadj_lowmem,&
36 ls_scf_submatrix_sign_ns
37 USE kinds, ONLY: dp,&
38 int_8
39 USE machine, ONLY: m_flush,&
40 m_walltime
41 USE mathconstants, ONLY: ifac
42 USE mathlib, ONLY: abnormal_value
43 USE message_passing, ONLY: mp_comm_type
44 USE submatrix_dissection, ONLY: submatrix_dissection_type
45#include "./base/base_uses.f90"
46
47 IMPLICIT NONE
48
49 PRIVATE
50
51 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'iterate_matrix'
52
53 TYPE :: eigbuf
54 REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: eigvals
55 REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: eigvecs
56 END TYPE eigbuf
57
58 INTERFACE purify_mcweeny
59 MODULE PROCEDURE purify_mcweeny_orth, purify_mcweeny_nonorth
60 END INTERFACE
61
62 PUBLIC :: invert_hotelling, matrix_sign_newton_schulz, matrix_sqrt_newton_schulz, &
63 matrix_sqrt_proot, matrix_sign_proot, matrix_sign_submatrix, matrix_exponential, &
64 matrix_sign_submatrix_mu_adjust, purify_mcweeny, invert_taylor, determinant
65
66CONTAINS
67
68! *****************************************************************************
69!> \brief Computes the determinant of a symmetric positive definite matrix
70!> using the trace of the matrix logarithm via Mercator series:
71!> det(A) = det(S)det(I+X)det(S), where S=diag(sqrt(Aii),..,sqrt(Ann))
72!> det(I+X) = Exp(Trace(Ln(I+X)))
73!> Ln(I+X) = X - X^2/2 + X^3/3 - X^4/4 + ..
74!> The series converges only if the Frobenius norm of X is less than 1.
75!> If it is more than one we compute (recursevily) the determinant of
76!> the square root of (I+X).
77!> \param matrix ...
78!> \param det - determinant
79!> \param threshold ...
80!> \par History
81!> 2015.04 created [Rustam Z Khaliullin]
82!> \author Rustam Z. Khaliullin
83! **************************************************************************************************
84 RECURSIVE SUBROUTINE determinant(matrix, det, threshold)
85
86 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
87 REAL(kind=dp), INTENT(INOUT) :: det
88 REAL(kind=dp), INTENT(IN) :: threshold
89
90 CHARACTER(LEN=*), PARAMETER :: routinen = 'determinant'
91
92 INTEGER :: handle, i, max_iter_lanczos, nsize, &
93 order_lanczos, sign_iter, unit_nr
94 INTEGER(KIND=int_8) :: flop1
95 INTEGER, SAVE :: recursion_depth = 0
96 REAL(kind=dp) :: det0, eps_lanczos, frobnorm, maxnorm, &
97 occ_matrix, t1, t2, trace
98 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: diagonal
99 TYPE(cp_logger_type), POINTER :: logger
100 TYPE(dbcsr_type) :: tmp1, tmp2, tmp3
101
102 CALL timeset(routinen, handle)
103
104 ! get a useful output_unit
105 logger => cp_get_default_logger()
106 IF (logger%para_env%is_source()) THEN
107 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
108 ELSE
109 unit_nr = -1
110 END IF
111
112 ! Note: tmp1 and tmp2 have the same matrix type as the
113 ! initial matrix (tmp3 does not have symmetry constraints)
114 ! this might lead to uninteded results with anti-symmetric
115 ! matrices
116 CALL dbcsr_create(tmp1, template=matrix, &
117 matrix_type=dbcsr_type_no_symmetry)
118 CALL dbcsr_create(tmp2, template=matrix, &
119 matrix_type=dbcsr_type_no_symmetry)
120 CALL dbcsr_create(tmp3, template=matrix, &
121 matrix_type=dbcsr_type_no_symmetry)
122
123 ! compute the product of the diagonal elements
124 block
125 TYPE(mp_comm_type) :: group
126 CALL dbcsr_get_info(matrix, nfullrows_total=nsize, group=group)
127 ALLOCATE (diagonal(nsize))
128 CALL dbcsr_get_diag(matrix, diagonal)
129 CALL group%sum(diagonal)
130 det = product(diagonal)
131 END block
132
133 ! create diagonal SQRTI matrix
134 diagonal(:) = 1.0_dp/(sqrt(diagonal(:)))
135 !ROLL CALL dbcsr_copy(tmp1,matrix)
136 CALL dbcsr_desymmetrize(matrix, tmp1)
137 CALL dbcsr_set(tmp1, 0.0_dp)
138 CALL dbcsr_set_diag(tmp1, diagonal)
139 CALL dbcsr_filter(tmp1, threshold)
140 DEALLOCATE (diagonal)
141
142 ! normalize the main diagonal, off-diagonal elements are scaled to
143 ! make the norm of the matrix less than 1
144 CALL dbcsr_multiply("N", "N", 1.0_dp, &
145 matrix, &
146 tmp1, &
147 0.0_dp, tmp3, &
148 filter_eps=threshold)
149 CALL dbcsr_multiply("N", "N", 1.0_dp, &
150 tmp1, &
151 tmp3, &
152 0.0_dp, tmp2, &
153 filter_eps=threshold)
154
155 ! subtract the main diagonal to create matrix X
156 CALL dbcsr_add_on_diag(tmp2, -1.0_dp)
157 frobnorm = dbcsr_frobenius_norm(tmp2)
158 IF (unit_nr > 0) THEN
159 IF (recursion_depth == 0) THEN
160 WRITE (unit_nr, '()')
161 ELSE
162 WRITE (unit_nr, '(T6,A28,1X,I15)') &
163 "Recursive iteration:", recursion_depth
164 END IF
165 WRITE (unit_nr, '(T6,A28,1X,F15.10)') &
166 "Frobenius norm:", frobnorm
167 CALL m_flush(unit_nr)
168 END IF
169
170 IF (frobnorm >= 1.0_dp) THEN
171
172 CALL dbcsr_add_on_diag(tmp2, 1.0_dp)
173 ! these controls should be provided as input
174 order_lanczos = 3
175 eps_lanczos = 1.0e-4_dp
176 max_iter_lanczos = 40
177 CALL matrix_sqrt_newton_schulz( &
178 tmp3, & ! output sqrt
179 tmp1, & ! output sqrti
180 tmp2, & ! input original
181 threshold=threshold, &
182 order=order_lanczos, &
183 eps_lanczos=eps_lanczos, &
184 max_iter_lanczos=max_iter_lanczos)
185 recursion_depth = recursion_depth + 1
186 CALL determinant(tmp3, det0, threshold)
187 recursion_depth = recursion_depth - 1
188 det = det*det0*det0
189
190 ELSE
191
192 ! create accumulator
193 CALL dbcsr_copy(tmp1, tmp2)
194 ! re-create to make use of symmetry
195 !ROLL CALL dbcsr_create(tmp3,template=matrix)
196
197 IF (unit_nr > 0) WRITE (unit_nr, *)
198
199 ! initialize the sign of the term
200 sign_iter = -1
201 DO i = 1, 100
202
203 t1 = m_walltime()
204
205 ! multiply X^i by X
206 ! note that the first iteration evaluates X^2
207 ! because the trace of X^1 is zero by construction
208 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp1, tmp2, &
209 0.0_dp, tmp3, &
210 filter_eps=threshold, &
211 flop=flop1)
212 CALL dbcsr_copy(tmp1, tmp3)
213
214 ! get trace
215 CALL dbcsr_trace(tmp1, trace)
216 trace = trace*sign_iter/(1.0_dp*(i + 1))
217 sign_iter = -sign_iter
218
219 ! update the determinant
220 det = det*exp(trace)
221
222 occ_matrix = dbcsr_get_occupation(tmp1)
223 maxnorm = dbcsr_maxabs(tmp1)
224
225 t2 = m_walltime()
226
227 IF (unit_nr > 0) THEN
228 WRITE (unit_nr, '(T6,A,1X,I3,1X,F7.5,F16.10,F10.3,F11.3)') &
229 "Determinant iter", i, occ_matrix, &
230 det, t2 - t1, &
231 flop1/(1.0e6_dp*max(0.001_dp, t2 - t1))
232 CALL m_flush(unit_nr)
233 END IF
234
235 ! exit if the trace is close to zero
236 IF (maxnorm < threshold) EXIT
237
238 END DO ! end iterations
239
240 IF (unit_nr > 0) THEN
241 WRITE (unit_nr, '()')
242 CALL m_flush(unit_nr)
243 END IF
244
245 END IF ! decide to do sqrt or not
246
247 IF (unit_nr > 0) THEN
248 IF (recursion_depth == 0) THEN
249 WRITE (unit_nr, '(T6,A28,1X,F15.10)') &
250 "Final determinant:", det
251 WRITE (unit_nr, '()')
252 ELSE
253 WRITE (unit_nr, '(T6,A28,1X,F15.10)') &
254 "Recursive determinant:", det
255 END IF
256 CALL m_flush(unit_nr)
257 END IF
258
259 CALL dbcsr_release(tmp1)
260 CALL dbcsr_release(tmp2)
261 CALL dbcsr_release(tmp3)
262
263 CALL timestop(handle)
264
265 END SUBROUTINE determinant
266
267! **************************************************************************************************
268!> \brief invert a symmetric positive definite diagonally dominant matrix
269!> \param matrix_inverse ...
270!> \param matrix ...
271!> \param threshold convergence threshold nased on the max abs
272!> \param use_inv_as_guess logical whether input can be used as guess for inverse
273!> \param norm_convergence convergence threshold for the 2-norm, useful for approximate solutions
274!> \param filter_eps filter_eps for matrix multiplications, if not passed nothing is filteres
275!> \param accelerator_order ...
276!> \param max_iter_lanczos ...
277!> \param eps_lanczos ...
278!> \param silent ...
279!> \par History
280!> 2010.10 created [Joost VandeVondele]
281!> 2011.10 guess option added [Rustam Z Khaliullin]
282!> \author Joost VandeVondele
283! **************************************************************************************************
284 SUBROUTINE invert_taylor(matrix_inverse, matrix, threshold, use_inv_as_guess, &
285 norm_convergence, filter_eps, accelerator_order, &
286 max_iter_lanczos, eps_lanczos, silent)
287
288 TYPE(dbcsr_type), INTENT(INOUT), TARGET :: matrix_inverse, matrix
289 REAL(kind=dp), INTENT(IN) :: threshold
290 LOGICAL, INTENT(IN), OPTIONAL :: use_inv_as_guess
291 REAL(kind=dp), INTENT(IN), OPTIONAL :: norm_convergence, filter_eps
292 INTEGER, INTENT(IN), OPTIONAL :: accelerator_order, max_iter_lanczos
293 REAL(kind=dp), INTENT(IN), OPTIONAL :: eps_lanczos
294 LOGICAL, INTENT(IN), OPTIONAL :: silent
295
296 CHARACTER(LEN=*), PARAMETER :: routinen = 'invert_Taylor'
297
298 INTEGER :: accelerator_type, handle, i, &
299 my_max_iter_lanczos, nrows, unit_nr
300 INTEGER(KIND=int_8) :: flop2
301 LOGICAL :: converged, use_inv_guess
302 REAL(kind=dp) :: coeff, convergence, maxnorm_matrix, &
303 my_eps_lanczos, occ_matrix, t1, t2
304 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: p_diagonal
305 TYPE(cp_logger_type), POINTER :: logger
306 TYPE(dbcsr_type), TARGET :: tmp1, tmp2, tmp3_sym
307
308 CALL timeset(routinen, handle)
309
310 logger => cp_get_default_logger()
311 IF (logger%para_env%is_source()) THEN
312 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
313 ELSE
314 unit_nr = -1
315 END IF
316 IF (PRESENT(silent)) THEN
317 IF (silent) unit_nr = -1
318 END IF
319
320 convergence = threshold
321 IF (PRESENT(norm_convergence)) convergence = norm_convergence
322
323 accelerator_type = 0
324 IF (PRESENT(accelerator_order)) accelerator_type = accelerator_order
325 IF (accelerator_type > 1) accelerator_type = 1
326
327 use_inv_guess = .false.
328 IF (PRESENT(use_inv_as_guess)) use_inv_guess = use_inv_as_guess
329
330 my_max_iter_lanczos = 64
331 my_eps_lanczos = 1.0e-3_dp
332 IF (PRESENT(max_iter_lanczos)) my_max_iter_lanczos = max_iter_lanczos
333 IF (PRESENT(eps_lanczos)) my_eps_lanczos = eps_lanczos
334
335 CALL dbcsr_create(tmp1, template=matrix_inverse, matrix_type=dbcsr_type_no_symmetry)
336 CALL dbcsr_create(tmp2, template=matrix_inverse, matrix_type=dbcsr_type_no_symmetry)
337 CALL dbcsr_create(tmp3_sym, template=matrix_inverse)
338
339 CALL dbcsr_get_info(matrix, nfullrows_total=nrows)
340 ALLOCATE (p_diagonal(nrows))
341
342 ! generate the initial guess
343 IF (.NOT. use_inv_guess) THEN
344
345 SELECT CASE (accelerator_type)
346 CASE (0)
347 ! use tmp1 to hold off-diagonal elements
348 CALL dbcsr_desymmetrize(matrix, tmp1)
349 p_diagonal(:) = 0.0_dp
350 CALL dbcsr_set_diag(tmp1, p_diagonal)
351 !CALL dbcsr_print(tmp1)
352 ! invert the main diagonal
353 CALL dbcsr_get_diag(matrix, p_diagonal)
354 DO i = 1, nrows
355 IF (p_diagonal(i) /= 0.0_dp) THEN
356 p_diagonal(i) = 1.0_dp/p_diagonal(i)
357 END IF
358 END DO
359 CALL dbcsr_set(matrix_inverse, 0.0_dp)
360 CALL dbcsr_add_on_diag(matrix_inverse, 1.0_dp)
361 CALL dbcsr_set_diag(matrix_inverse, p_diagonal)
362 CASE DEFAULT
363 cpabort("Illegal accelerator order")
364 END SELECT
365
366 ELSE
367
368 cpabort("Guess is NYI")
369
370 END IF
371
372 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp1, matrix_inverse, &
373 0.0_dp, tmp2, filter_eps=filter_eps)
374
375 IF (unit_nr > 0) WRITE (unit_nr, *)
376
377 ! scale the approximate inverse to be within the convergence radius
378 t1 = m_walltime()
379
380 ! done with the initial guess, start iterations
381 converged = .false.
382 CALL dbcsr_desymmetrize(matrix_inverse, tmp1)
383 coeff = 1.0_dp
384 DO i = 1, 100
385
386 ! coeff = +/- 1
387 coeff = -1.0_dp*coeff
388 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp1, tmp2, 0.0_dp, &
389 tmp3_sym, &
390 flop=flop2, filter_eps=filter_eps)
391 !flop=flop2)
392 CALL dbcsr_add(matrix_inverse, tmp3_sym, 1.0_dp, coeff)
393 CALL dbcsr_release(tmp1)
394 CALL dbcsr_create(tmp1, template=matrix_inverse, matrix_type=dbcsr_type_no_symmetry)
395 CALL dbcsr_desymmetrize(tmp3_sym, tmp1)
396
397 ! for the convergence check
398 maxnorm_matrix = dbcsr_maxabs(tmp3_sym)
399
400 t2 = m_walltime()
401 occ_matrix = dbcsr_get_occupation(matrix_inverse)
402
403 IF (unit_nr > 0) THEN
404 WRITE (unit_nr, '(T6,A,1X,I3,1X,F10.8,E12.3,F12.3,F13.3)') "Taylor iter", i, occ_matrix, &
405 maxnorm_matrix, t2 - t1, &
406 flop2/(1.0e6_dp*max(0.001_dp, t2 - t1))
407 CALL m_flush(unit_nr)
408 END IF
409
410 IF (maxnorm_matrix < convergence) THEN
411 converged = .true.
412 EXIT
413 END IF
414
415 t1 = m_walltime()
416
417 END DO
418
419 !last convergence check
420 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, matrix_inverse, 0.0_dp, tmp1, &
421 filter_eps=filter_eps)
422 CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
423 !frob_matrix = dbcsr_frobenius_norm(tmp1)
424 maxnorm_matrix = dbcsr_maxabs(tmp1)
425 IF (unit_nr > 0) THEN
426 WRITE (unit_nr, '(T6,A,E12.5)') "Final Taylor error", maxnorm_matrix
427 WRITE (unit_nr, '()')
428 CALL m_flush(unit_nr)
429 END IF
430 IF (maxnorm_matrix > convergence) THEN
431 converged = .false.
432 IF (unit_nr > 0) THEN
433 WRITE (unit_nr, *) 'Final convergence check failed'
434 END IF
435 END IF
436
437 IF (.NOT. converged) THEN
438 cpabort("Taylor inversion did not converge")
439 END IF
440
441 CALL dbcsr_release(tmp1)
442 CALL dbcsr_release(tmp2)
443 CALL dbcsr_release(tmp3_sym)
444
445 DEALLOCATE (p_diagonal)
446
447 CALL timestop(handle)
448
449 END SUBROUTINE invert_taylor
450
451! **************************************************************************************************
452!> \brief invert a symmetric positive definite matrix by Hotelling's method
453!> explicit symmetrization makes this code not suitable for other matrix types
454!> Currently a bit messy with the options, to to be cleaned soon
455!> \param matrix_inverse ...
456!> \param matrix ...
457!> \param threshold convergence threshold nased on the max abs
458!> \param use_inv_as_guess logical whether input can be used as guess for inverse
459!> \param norm_convergence convergence threshold for the 2-norm, useful for approximate solutions
460!> \param filter_eps filter_eps for matrix multiplications, if not passed nothing is filteres
461!> \param accelerator_order ...
462!> \param max_iter_lanczos ...
463!> \param eps_lanczos ...
464!> \param silent ...
465!> \par History
466!> 2010.10 created [Joost VandeVondele]
467!> 2011.10 guess option added [Rustam Z Khaliullin]
468!> \author Joost VandeVondele
469! **************************************************************************************************
470 SUBROUTINE invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, &
471 norm_convergence, filter_eps, accelerator_order, &
472 max_iter_lanczos, eps_lanczos, silent)
473
474 TYPE(dbcsr_type), INTENT(INOUT), TARGET :: matrix_inverse, matrix
475 REAL(kind=dp), INTENT(IN) :: threshold
476 LOGICAL, INTENT(IN), OPTIONAL :: use_inv_as_guess
477 REAL(kind=dp), INTENT(IN), OPTIONAL :: norm_convergence, filter_eps
478 INTEGER, INTENT(IN), OPTIONAL :: accelerator_order, max_iter_lanczos
479 REAL(kind=dp), INTENT(IN), OPTIONAL :: eps_lanczos
480 LOGICAL, INTENT(IN), OPTIONAL :: silent
481
482 CHARACTER(LEN=*), PARAMETER :: routinen = 'invert_Hotelling'
483
484 INTEGER :: accelerator_type, handle, i, &
485 my_max_iter_lanczos, unit_nr
486 INTEGER(KIND=int_8) :: flop1, flop2
487 LOGICAL :: arnoldi_converged, converged, &
488 use_inv_guess
489 REAL(kind=dp) :: convergence, frob_matrix, gershgorin_norm, max_ev, maxnorm_matrix, min_ev, &
490 my_eps_lanczos, my_filter_eps, occ_matrix, scalingf, t1, t2
491 TYPE(cp_logger_type), POINTER :: logger
492 TYPE(dbcsr_type), TARGET :: tmp1, tmp2
493
494 !TYPE(arnoldi_env_type) :: arnoldi_env
495 !TYPE(dbcsr_p_type), DIMENSION(1) :: mymat
496
497 CALL timeset(routinen, handle)
498
499 logger => cp_get_default_logger()
500 IF (logger%para_env%is_source()) THEN
501 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
502 ELSE
503 unit_nr = -1
504 END IF
505 IF (PRESENT(silent)) THEN
506 IF (silent) unit_nr = -1
507 END IF
508
509 convergence = threshold
510 IF (PRESENT(norm_convergence)) convergence = norm_convergence
511
512 accelerator_type = 1
513 IF (PRESENT(accelerator_order)) accelerator_type = accelerator_order
514 IF (accelerator_type > 1) accelerator_type = 1
515
516 use_inv_guess = .false.
517 IF (PRESENT(use_inv_as_guess)) use_inv_guess = use_inv_as_guess
518
519 my_max_iter_lanczos = 64
520 my_eps_lanczos = 1.0e-3_dp
521 IF (PRESENT(max_iter_lanczos)) my_max_iter_lanczos = max_iter_lanczos
522 IF (PRESENT(eps_lanczos)) my_eps_lanczos = eps_lanczos
523
524 my_filter_eps = threshold
525 IF (PRESENT(filter_eps)) my_filter_eps = filter_eps
526
527 ! generate the initial guess
528 IF (.NOT. use_inv_guess) THEN
529
530 SELECT CASE (accelerator_type)
531 CASE (0)
532 gershgorin_norm = dbcsr_gershgorin_norm(matrix)
533 frob_matrix = dbcsr_frobenius_norm(matrix)
534 CALL dbcsr_set(matrix_inverse, 0.0_dp)
535 CALL dbcsr_add_on_diag(matrix_inverse, 1/min(gershgorin_norm, frob_matrix))
536 CASE (1)
537 ! initialize matrix to unity and use arnoldi (below) to scale it into the convergence range
538 CALL dbcsr_set(matrix_inverse, 0.0_dp)
539 CALL dbcsr_add_on_diag(matrix_inverse, 1.0_dp)
540 CASE DEFAULT
541 cpabort("Illegal accelerator order")
542 END SELECT
543
544 ! everything commutes, therefore our all products will be symmetric
545 CALL dbcsr_create(tmp1, template=matrix_inverse)
546
547 ELSE
548
549 ! It is unlikely that our guess will commute with the matrix, therefore the first product will
550 ! be non symmetric
551 CALL dbcsr_create(tmp1, template=matrix_inverse, matrix_type=dbcsr_type_no_symmetry)
552
553 END IF
554
555 CALL dbcsr_create(tmp2, template=matrix_inverse)
556
557 IF (unit_nr > 0) WRITE (unit_nr, *)
558
559 ! scale the approximate inverse to be within the convergence radius
560 t1 = m_walltime()
561
562 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_inverse, matrix, &
563 0.0_dp, tmp1, flop=flop1, filter_eps=my_filter_eps)
564
565 IF (accelerator_type == 1) THEN
566
567 ! scale the matrix to get into the convergence range
568 CALL arnoldi_extremal(tmp1, max_ev, min_ev, threshold=my_eps_lanczos, &
569 max_iter=my_max_iter_lanczos, converged=arnoldi_converged)
570 !mymat(1)%matrix => tmp1
571 !CALL setup_arnoldi_env(arnoldi_env, mymat, max_iter=30, threshold=1.0E-3_dp, selection_crit=1, &
572 ! nval_request=2, nrestarts=2, generalized_ev=.FALSE., iram=.TRUE.)
573 !CALL arnoldi_ev(mymat, arnoldi_env)
574 !max_eV = REAL(get_selected_ritz_val(arnoldi_env, 2), dp)
575 !min_eV = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
576 !CALL deallocate_arnoldi_env(arnoldi_env)
577
578 IF (unit_nr > 0) THEN
579 WRITE (unit_nr, *)
580 WRITE (unit_nr, '(T6,A,1X,L1,A,E12.3)') "Lanczos converged: ", arnoldi_converged, " threshold:", my_eps_lanczos
581 WRITE (unit_nr, '(T6,A,1X,E12.3,E12.3)') "Est. extremal eigenvalues:", max_ev, min_ev
582 WRITE (unit_nr, '(T6,A,1X,E12.3)') "Est. condition number :", max_ev/max(min_ev, epsilon(min_ev))
583 END IF
584
585 ! 2.0 would be the correct scaling however, we should make sure here, that we are in the convergence radius
586 scalingf = 1.9_dp/(max_ev + min_ev)
587 CALL dbcsr_scale(tmp1, scalingf)
588 CALL dbcsr_scale(matrix_inverse, scalingf)
589 min_ev = min_ev*scalingf
590
591 END IF
592
593 ! done with the initial guess, start iterations
594 converged = .false.
595 DO i = 1, 100
596
597 ! tmp1 = S^-1 S
598
599 ! for the convergence check
600 CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
601 maxnorm_matrix = dbcsr_maxabs(tmp1)
602 CALL dbcsr_add_on_diag(tmp1, +1.0_dp)
603
604 ! tmp2 = S^-1 S S^-1
605 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp1, matrix_inverse, 0.0_dp, tmp2, &
606 flop=flop2, filter_eps=my_filter_eps)
607 ! S^-1_{n+1} = 2 S^-1 - S^-1 S S^-1
608 CALL dbcsr_add(matrix_inverse, tmp2, 2.0_dp, -1.0_dp)
609
610 CALL dbcsr_filter(matrix_inverse, my_filter_eps)
611 t2 = m_walltime()
612 occ_matrix = dbcsr_get_occupation(matrix_inverse)
613
614 ! use the scalar form of the algorithm to trace the EV
615 IF (accelerator_type == 1) THEN
616 min_ev = min_ev*(2.0_dp - min_ev)
617 IF (PRESENT(norm_convergence)) maxnorm_matrix = abs(min_ev - 1.0_dp)
618 END IF
619
620 IF (unit_nr > 0) THEN
621 WRITE (unit_nr, '(T6,A,1X,I3,1X,F10.8,E12.3,F12.3,F13.3)') "Hotelling iter", i, occ_matrix, &
622 maxnorm_matrix, t2 - t1, &
623 (flop1 + flop2)/(1.0e6_dp*max(0.001_dp, t2 - t1))
624 CALL m_flush(unit_nr)
625 END IF
626
627 IF (maxnorm_matrix < convergence) THEN
628 converged = .true.
629 EXIT
630 END IF
631
632 ! scale the matrix for improved convergence
633 IF (accelerator_type == 1) THEN
634 min_ev = min_ev*2.0_dp/(min_ev + 1.0_dp)
635 CALL dbcsr_scale(matrix_inverse, 2.0_dp/(min_ev + 1.0_dp))
636 END IF
637
638 t1 = m_walltime()
639 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_inverse, matrix, &
640 0.0_dp, tmp1, flop=flop1, filter_eps=my_filter_eps)
641
642 END DO
643
644 IF (.NOT. converged) THEN
645 cpabort("Hotelling inversion did not converge")
646 END IF
647
648 ! try to symmetrize the output matrix
649 IF (dbcsr_get_matrix_type(matrix_inverse) == dbcsr_type_no_symmetry) THEN
650 CALL dbcsr_transposed(tmp2, matrix_inverse)
651 CALL dbcsr_add(matrix_inverse, tmp2, 0.5_dp, 0.5_dp)
652 END IF
653
654 IF (unit_nr > 0) THEN
655! WRITE(unit_nr,'(T6,A,1X,I3,1X,F10.8,E12.3)') "Final Hotelling ",i,occ_matrix,&
656! !frob_matrix/frob_matrix_base
657! maxnorm_matrix
658 WRITE (unit_nr, '()')
659 CALL m_flush(unit_nr)
660 END IF
661
662 CALL dbcsr_release(tmp1)
663 CALL dbcsr_release(tmp2)
664
665 CALL timestop(handle)
666
667 END SUBROUTINE invert_hotelling
668
669! **************************************************************************************************
670!> \brief compute the sign a matrix using Newton-Schulz iterations
671!> \param matrix_sign ...
672!> \param matrix ...
673!> \param threshold ...
674!> \param sign_order ...
675!> \param iounit ...
676!> \par History
677!> 2010.10 created [Joost VandeVondele]
678!> 2019.05 extended to order byxond 2 [Robert Schade]
679!> \author Joost VandeVondele, Robert Schade
680! **************************************************************************************************
681 SUBROUTINE matrix_sign_newton_schulz(matrix_sign, matrix, threshold, sign_order, iounit)
682
683 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_sign, matrix
684 REAL(kind=dp), INTENT(IN) :: threshold
685 INTEGER, INTENT(IN), OPTIONAL :: sign_order, iounit
686
687 CHARACTER(LEN=*), PARAMETER :: routinen = 'matrix_sign_Newton_Schulz'
688
689 INTEGER :: count, handle, i, order, unit_nr
690 INTEGER(KIND=int_8) :: flops
691 REAL(kind=dp) :: a0, a1, a2, a3, a4, a5, floptot, &
692 frob_matrix, frob_matrix_base, &
693 gersh_matrix, occ_matrix, prefactor, &
694 t1, t2
695 TYPE(cp_logger_type), POINTER :: logger
696 TYPE(dbcsr_type) :: tmp1, tmp2, tmp3, tmp4
697
698 CALL timeset(routinen, handle)
699
700 IF (PRESENT(iounit)) THEN
701 unit_nr = iounit
702 ELSE
703 logger => cp_get_default_logger()
704 IF (logger%para_env%is_source()) THEN
705 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
706 ELSE
707 unit_nr = -1
708 END IF
709 END IF
710
711 IF (PRESENT(sign_order)) THEN
712 order = sign_order
713 ELSE
714 order = 2
715 END IF
716
717 CALL dbcsr_create(tmp1, template=matrix_sign)
718
719 CALL dbcsr_create(tmp2, template=matrix_sign)
720 IF (abs(order) >= 4) THEN
721 CALL dbcsr_create(tmp3, template=matrix_sign)
722 END IF
723 IF (abs(order) > 4) THEN
724 CALL dbcsr_create(tmp4, template=matrix_sign)
725 END IF
726
727 CALL dbcsr_copy(matrix_sign, matrix)
728 CALL dbcsr_filter(matrix_sign, threshold)
729
730 ! scale the matrix to get into the convergence range
731 frob_matrix = dbcsr_frobenius_norm(matrix_sign)
732 gersh_matrix = dbcsr_gershgorin_norm(matrix_sign)
733 CALL dbcsr_scale(matrix_sign, 1/min(frob_matrix, gersh_matrix))
734
735 IF (unit_nr > 0) WRITE (unit_nr, *)
736
737 count = 0
738 DO i = 1, 100
739 floptot = 0_dp
740 t1 = m_walltime()
741 ! tmp1 = X * X
742 CALL dbcsr_multiply("N", "N", -1.0_dp, matrix_sign, matrix_sign, 0.0_dp, tmp1, &
743 filter_eps=threshold, flop=flops)
744 floptot = floptot + flops
745
746 ! check convergence (frob norm of what should be the identity matrix minus identity matrix)
747 frob_matrix_base = dbcsr_frobenius_norm(tmp1)
748 CALL dbcsr_add_on_diag(tmp1, +1.0_dp)
749 frob_matrix = dbcsr_frobenius_norm(tmp1)
750
751 ! f(y) approx 1/sqrt(1-y)
752 ! f(y)=1+y/2+3/8*y^2+5/16*y^3+35/128*y^4+63/256*y^5+231/1024*y^6
753 ! f2(y)=1+y/2=1/2*(2+y)
754 ! f3(y)=1+y/2+3/8*y^2=3/8*(8/3+4/3*y+y^2)
755 ! f4(y)=1+y/2+3/8*y^2+5/16*y^3=5/16*(16/5+8/5*y+6/5*y^2+y^3)
756 ! f5(y)=1+y/2+3/8*y^2+5/16*y^3+35/128*y^4=35/128*(128/35+128/70*y+48/35*y^2+8/7*y^3+y^4)
757 ! z(y)=(y+a_0)*y+a_1
758 ! f5(y)=35/128*((z(y)+y+a_2)*z(y)+a_3)
759 ! =35/128*((a_1^2+a_1a_2+a_3)+(2*a_0a_1+a_1+a_0a_2)y+(a_0^2+a_0+2a_1+a_2)y^2+(2a_0+1)y^3+y^4)
760 ! a_0=1/14
761 ! a_1=23819/13720
762 ! a_2=1269/980-2a_1=-3734/1715
763 ! a_3=832591127/188238400
764 ! f6(y)=1+y/2+3/8*y^2+5/16*y^3+35/128*y^4+63/256*y^5
765 ! =63/256*(256/63 + (128 y)/63 + (32 y^2)/21 + (80 y^3)/63 + (10 y^4)/9 + y^5)
766 ! f7(y)=1+y/2+3/8*y^2+5/16*y^3+35/128*y^4+63/256*y^5+231/1024*y^6
767 ! =231/1024*(1024/231+512/231*y+128/77*y^2+320/231*y^3+40/33*y^4+12/11*y^5+y^6)
768 ! z(y)=(y+a_0)*y+a_1
769 ! w(y)=(y+a_2)*z(y)+a_3
770 ! f7(y)=(w(y)+z(y)+a_4)*w(y)+a_5
771 ! a_0= 1.3686502058092053653287666647611728507211996691324048468010382350359929055186612505791532871573242422
772 ! a_1= 1.7089671854477436685850554669524985556296280184497503489303331821456795715195510972774979091893741568
773 ! a_2=-1.3231956603546599107833121193066273961757451236778593922555836895814474509732067051246078326118696968
774 ! a_3= 3.9876642330847931291749479958277754186675336169578593000744380254770411483327581042259415937710270453
775 ! a_4=-3.7273299006476825027065704937541279833880400042556351139273912137942678919776364526511485025132991667
776 ! a_5= 4.9369932474103023792021351907971943220607580694533770325967170245194362399287150565595441897740173578
777 !
778 ! y=1-X*X
779
780 ! tmp1 = I-x*x
781 IF (order == 2) THEN
782 prefactor = 0.5_dp
783
784 ! update the above to 3*I-X*X
785 CALL dbcsr_add_on_diag(tmp1, +2.0_dp)
786 occ_matrix = dbcsr_get_occupation(matrix_sign)
787 ELSE IF (order == 3) THEN
788 ! with one multiplication
789 ! tmp1=y
790 CALL dbcsr_copy(tmp2, tmp1)
791 CALL dbcsr_scale(tmp1, 4.0_dp/3.0_dp)
792 CALL dbcsr_add_on_diag(tmp1, 8.0_dp/3.0_dp)
793
794 ! tmp2=y^2
795 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp2, tmp2, 1.0_dp, tmp1, &
796 filter_eps=threshold, flop=flops)
797 floptot = floptot + flops
798 prefactor = 3.0_dp/8.0_dp
799
800 ELSE IF (order == 4) THEN
801 ! with two multiplications
802 ! tmp1=y
803 CALL dbcsr_copy(tmp3, tmp1)
804 CALL dbcsr_scale(tmp1, 8.0_dp/5.0_dp)
805 CALL dbcsr_add_on_diag(tmp1, 16.0_dp/5.0_dp)
806
807 !
808 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp3, tmp3, 0.0_dp, tmp2, &
809 filter_eps=threshold, flop=flops)
810 floptot = floptot + flops
811
812 CALL dbcsr_add(tmp1, tmp2, 1.0_dp, 6.0_dp/5.0_dp)
813
814 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp2, tmp3, 1.0_dp, tmp1, &
815 filter_eps=threshold, flop=flops)
816 floptot = floptot + flops
817
818 prefactor = 5.0_dp/16.0_dp
819 ELSE IF (order == -5) THEN
820 ! with three multiplications
821 ! tmp1=y
822 CALL dbcsr_copy(tmp3, tmp1)
823 CALL dbcsr_scale(tmp1, 128.0_dp/70.0_dp)
824 CALL dbcsr_add_on_diag(tmp1, 128.0_dp/35.0_dp)
825
826 !
827 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp3, tmp3, 0.0_dp, tmp2, &
828 filter_eps=threshold, flop=flops)
829 floptot = floptot + flops
830
831 CALL dbcsr_add(tmp1, tmp2, 1.0_dp, 48.0_dp/35.0_dp)
832
833 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp2, tmp3, 0.0_dp, tmp4, &
834 filter_eps=threshold, flop=flops)
835 floptot = floptot + flops
836
837 CALL dbcsr_add(tmp1, tmp4, 1.0_dp, 8.0_dp/7.0_dp)
838
839 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp4, tmp3, 1.0_dp, tmp1, &
840 filter_eps=threshold, flop=flops)
841 floptot = floptot + flops
842
843 prefactor = 35.0_dp/128.0_dp
844 ELSE IF (order == 5) THEN
845 ! with two multiplications
846 ! z(y)=(y+a_0)*y+a_1
847 ! f5(y)=35/128*((z(y)+y+a_2)*z(y)+a_3)
848 ! =35/128*((a_1^2+a_1a_2+a_3)+(2*a_0a_1+a_1+a_0a_2)y+(a_0^2+a_0+2a_1+a_2)y^2+(2a_0+1)y^3+y^4)
849 ! a_0=1/14
850 ! a_1=23819/13720
851 ! a_2=1269/980-2a_1=-3734/1715
852 ! a_3=832591127/188238400
853 a0 = 1.0_dp/14.0_dp
854 a1 = 23819.0_dp/13720.0_dp
855 a2 = -3734_dp/1715.0_dp
856 a3 = 832591127_dp/188238400.0_dp
857
858 ! tmp1=y
859 ! tmp3=z
860 CALL dbcsr_copy(tmp3, tmp1)
861 CALL dbcsr_add_on_diag(tmp3, a0)
862 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp3, tmp1, 0.0_dp, tmp2, &
863 filter_eps=threshold, flop=flops)
864 floptot = floptot + flops
865 CALL dbcsr_add_on_diag(tmp2, a1)
866
867 CALL dbcsr_add_on_diag(tmp1, a2)
868 CALL dbcsr_add(tmp1, tmp2, 1.0_dp, 1.0_dp)
869 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp1, tmp2, 0.0_dp, tmp3, &
870 filter_eps=threshold, flop=flops)
871 floptot = floptot + flops
872 CALL dbcsr_add_on_diag(tmp3, a3)
873 CALL dbcsr_copy(tmp1, tmp3)
874
875 prefactor = 35.0_dp/128.0_dp
876 ELSE IF (order == 6) THEN
877 ! with four multiplications
878 ! f6(y)=63/256*(256/63 + (128 y)/63 + (32 y^2)/21 + (80 y^3)/63 + (10 y^4)/9 + y^5)
879 ! tmp1=y
880 CALL dbcsr_copy(tmp3, tmp1)
881 CALL dbcsr_scale(tmp1, 128.0_dp/63.0_dp)
882 CALL dbcsr_add_on_diag(tmp1, 256.0_dp/63.0_dp)
883
884 !
885 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp3, tmp3, 0.0_dp, tmp2, &
886 filter_eps=threshold, flop=flops)
887 floptot = floptot + flops
888
889 CALL dbcsr_add(tmp1, tmp2, 1.0_dp, 32.0_dp/21.0_dp)
890
891 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp2, tmp3, 0.0_dp, tmp4, &
892 filter_eps=threshold, flop=flops)
893 floptot = floptot + flops
894
895 CALL dbcsr_add(tmp1, tmp4, 1.0_dp, 80.0_dp/63.0_dp)
896
897 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp4, tmp3, 0.0_dp, tmp2, &
898 filter_eps=threshold, flop=flops)
899 floptot = floptot + flops
900
901 CALL dbcsr_add(tmp1, tmp2, 1.0_dp, 10.0_dp/9.0_dp)
902
903 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp2, tmp3, 1.0_dp, tmp1, &
904 filter_eps=threshold, flop=flops)
905 floptot = floptot + flops
906
907 prefactor = 63.0_dp/256.0_dp
908 ELSE IF (order == 7) THEN
909 ! with three multiplications
910
911 a0 = 1.3686502058092053653287666647611728507211996691324048468010382350359929055186612505791532871573242422_dp
912 a1 = 1.7089671854477436685850554669524985556296280184497503489303331821456795715195510972774979091893741568_dp
913 a2 = -1.3231956603546599107833121193066273961757451236778593922555836895814474509732067051246078326118696968_dp
914 a3 = 3.9876642330847931291749479958277754186675336169578593000744380254770411483327581042259415937710270453_dp
915 a4 = -3.7273299006476825027065704937541279833880400042556351139273912137942678919776364526511485025132991667_dp
916 a5 = 4.9369932474103023792021351907971943220607580694533770325967170245194362399287150565595441897740173578_dp
917 ! =231/1024*(1024/231+512/231*y+128/77*y^2+320/231*y^3+40/33*y^4+12/11*y^5+y^6)
918 ! z(y)=(y+a_0)*y+a_1
919 ! w(y)=(y+a_2)*z(y)+a_3
920 ! f7(y)=(w(y)+z(y)+a_4)*w(y)+a_5
921
922 ! tmp1=y
923 ! tmp3=z
924 CALL dbcsr_copy(tmp3, tmp1)
925 CALL dbcsr_add_on_diag(tmp3, a0)
926 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp3, tmp1, 0.0_dp, tmp2, &
927 filter_eps=threshold, flop=flops)
928 floptot = floptot + flops
929 CALL dbcsr_add_on_diag(tmp2, a1)
930
931 ! tmp4=w
932 CALL dbcsr_copy(tmp4, tmp1)
933 CALL dbcsr_add_on_diag(tmp4, a2)
934 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp4, tmp2, 0.0_dp, tmp3, &
935 filter_eps=threshold, flop=flops)
936 floptot = floptot + flops
937 CALL dbcsr_add_on_diag(tmp3, a3)
938
939 CALL dbcsr_add(tmp2, tmp3, 1.0_dp, 1.0_dp)
940 CALL dbcsr_add_on_diag(tmp2, a4)
941 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp2, tmp3, 0.0_dp, tmp1, &
942 filter_eps=threshold, flop=flops)
943 floptot = floptot + flops
944 CALL dbcsr_add_on_diag(tmp1, a5)
945
946 prefactor = 231.0_dp/1024.0_dp
947 ELSE
948 cpabort("requested order is not implemented.")
949 END IF
950
951 ! tmp2 = X * prefactor *
952 CALL dbcsr_multiply("N", "N", prefactor, matrix_sign, tmp1, 0.0_dp, tmp2, &
953 filter_eps=threshold, flop=flops)
954 floptot = floptot + flops
955
956 ! done iterating
957 ! CALL dbcsr_filter(tmp2,threshold)
958 CALL dbcsr_copy(matrix_sign, tmp2)
959 t2 = m_walltime()
960
961 occ_matrix = dbcsr_get_occupation(matrix_sign)
962
963 IF (unit_nr > 0) THEN
964 WRITE (unit_nr, '(T6,A,1X,I3,1X,F10.8,E12.3,F12.3,F13.3)') "NS sign iter ", i, occ_matrix, &
965 frob_matrix/frob_matrix_base, t2 - t1, &
966 floptot/(1.0e6_dp*max(0.001_dp, t2 - t1))
967 CALL m_flush(unit_nr)
968 END IF
969
970 ! frob_matrix/frob_matrix_base < SQRT(threshold)
971 IF (frob_matrix*frob_matrix < (threshold*frob_matrix_base*frob_matrix_base)) EXIT
972
973 END DO
974
975 ! this check is not really needed
976 CALL dbcsr_multiply("N", "N", +1.0_dp, matrix_sign, matrix_sign, 0.0_dp, tmp1, &
977 filter_eps=threshold)
978 frob_matrix_base = dbcsr_frobenius_norm(tmp1)
979 CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
980 frob_matrix = dbcsr_frobenius_norm(tmp1)
981 occ_matrix = dbcsr_get_occupation(matrix_sign)
982 IF (unit_nr > 0) THEN
983 WRITE (unit_nr, '(T6,A,1X,I3,1X,F10.8,E12.3)') "Final NS sign iter", i, occ_matrix, &
984 frob_matrix/frob_matrix_base
985 WRITE (unit_nr, '()')
986 CALL m_flush(unit_nr)
987 END IF
988
989 CALL dbcsr_release(tmp1)
990 CALL dbcsr_release(tmp2)
991 IF (abs(order) >= 4) THEN
992 CALL dbcsr_release(tmp3)
993 END IF
994 IF (abs(order) > 4) THEN
995 CALL dbcsr_release(tmp4)
996 END IF
997
998 CALL timestop(handle)
999
1000 END SUBROUTINE matrix_sign_newton_schulz
1001
1002 ! **************************************************************************************************
1003!> \brief compute the sign a matrix using the general algorithm for the p-th root of Richters et al.
1004!> Commun. Comput. Phys., 25 (2019), pp. 564-585.
1005!> \param matrix_sign ...
1006!> \param matrix ...
1007!> \param threshold ...
1008!> \param sign_order ...
1009!> \par History
1010!> 2019.03 created [Robert Schade]
1011!> \author Robert Schade
1012! **************************************************************************************************
1013 SUBROUTINE matrix_sign_proot(matrix_sign, matrix, threshold, sign_order)
1014
1015 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_sign, matrix
1016 REAL(kind=dp), INTENT(IN) :: threshold
1017 INTEGER, INTENT(IN), OPTIONAL :: sign_order
1018
1019 CHARACTER(LEN=*), PARAMETER :: routinen = 'matrix_sign_proot'
1020
1021 INTEGER :: handle, order, unit_nr
1022 INTEGER(KIND=int_8) :: flop0, flop1, flop2
1023 LOGICAL :: converged, symmetrize
1024 REAL(kind=dp) :: frob_matrix, frob_matrix_base, occ_matrix
1025 TYPE(cp_logger_type), POINTER :: logger
1026 TYPE(dbcsr_type) :: matrix2, matrix_sqrt, matrix_sqrt_inv, &
1027 tmp1, tmp2
1028
1029 CALL cite_reference(richters2018)
1030
1031 CALL timeset(routinen, handle)
1032
1033 logger => cp_get_default_logger()
1034 IF (logger%para_env%is_source()) THEN
1035 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1036 ELSE
1037 unit_nr = -1
1038 END IF
1039
1040 IF (PRESENT(sign_order)) THEN
1041 order = sign_order
1042 ELSE
1043 order = 2
1044 END IF
1045
1046 CALL dbcsr_create(tmp1, template=matrix_sign)
1047
1048 CALL dbcsr_create(tmp2, template=matrix_sign)
1049
1050 CALL dbcsr_create(matrix2, template=matrix, matrix_type=dbcsr_type_no_symmetry)
1051 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, matrix, 0.0_dp, matrix2, &
1052 filter_eps=threshold, flop=flop0)
1053 !CALL dbcsr_filter(matrix2, threshold)
1054
1055 !CALL dbcsr_copy(matrix_sign, matrix)
1056 !CALL dbcsr_filter(matrix_sign, threshold)
1057
1058 IF (unit_nr > 0) WRITE (unit_nr, *)
1059
1060 CALL dbcsr_create(matrix_sqrt, template=matrix2)
1061 CALL dbcsr_create(matrix_sqrt_inv, template=matrix2)
1062 IF (unit_nr > 0) WRITE (unit_nr, *) "Threshold=", threshold
1063
1064 symmetrize = .false.
1065 CALL matrix_sqrt_proot(matrix_sqrt, matrix_sqrt_inv, matrix2, threshold, order, &
1066 0.01_dp, 100, symmetrize, converged)
1067
1068 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, matrix_sqrt_inv, 0.0_dp, matrix_sign, &
1069 filter_eps=threshold, flop=flop1)
1070
1071 ! this check is not really needed
1072 CALL dbcsr_multiply("N", "N", +1.0_dp, matrix_sign, matrix_sign, 0.0_dp, tmp1, &
1073 filter_eps=threshold, flop=flop2)
1074 frob_matrix_base = dbcsr_frobenius_norm(tmp1)
1075 CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
1076 frob_matrix = dbcsr_frobenius_norm(tmp1)
1077 occ_matrix = dbcsr_get_occupation(matrix_sign)
1078 IF (unit_nr > 0) THEN
1079 WRITE (unit_nr, '(T6,A,F10.8,E12.3)') "Final proot sign iter", occ_matrix, &
1080 frob_matrix/frob_matrix_base
1081 WRITE (unit_nr, '()')
1082 CALL m_flush(unit_nr)
1083 END IF
1084
1085 CALL dbcsr_release(tmp1)
1086 CALL dbcsr_release(tmp2)
1087 CALL dbcsr_release(matrix2)
1088 CALL dbcsr_release(matrix_sqrt)
1089 CALL dbcsr_release(matrix_sqrt_inv)
1090
1091 CALL timestop(handle)
1092
1093 END SUBROUTINE matrix_sign_proot
1094
1095! **************************************************************************************************
1096!> \brief compute the sign of a dense matrix using Newton-Schulz iterations
1097!> \param matrix_sign ...
1098!> \param matrix ...
1099!> \param matrix_id ...
1100!> \param threshold ...
1101!> \param sign_order ...
1102!> \author Michael Lass, Robert Schade
1103! **************************************************************************************************
1104 SUBROUTINE dense_matrix_sign_newton_schulz(matrix_sign, matrix, matrix_id, threshold, sign_order)
1105
1106 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix_sign
1107 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: matrix
1108 INTEGER, INTENT(IN), OPTIONAL :: matrix_id
1109 REAL(kind=dp), INTENT(IN) :: threshold
1110 INTEGER, INTENT(IN), OPTIONAL :: sign_order
1111
1112 CHARACTER(LEN=*), PARAMETER :: routinen = 'dense_matrix_sign_Newton_Schulz'
1113
1114 INTEGER :: handle, i, j, sz, unit_nr
1115 LOGICAL :: converged
1116 REAL(kind=dp) :: frob_matrix, frob_matrix_base, &
1117 gersh_matrix, prefactor, scaling_factor
1118 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp1, tmp2
1119 REAL(kind=dp), DIMENSION(1) :: work
1120 REAL(kind=dp), EXTERNAL :: dlange
1121 TYPE(cp_logger_type), POINTER :: logger
1122
1123 CALL timeset(routinen, handle)
1124
1125 ! print output on all ranks
1126 logger => cp_get_default_logger()
1127 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1128
1129 ! scale the matrix to get into the convergence range
1130 sz = SIZE(matrix, 1)
1131 frob_matrix = dlange('F', sz, sz, matrix, sz, work) !dbcsr_frobenius_norm(matrix_sign)
1132 gersh_matrix = dlange('1', sz, sz, matrix, sz, work) !dbcsr_gershgorin_norm(matrix_sign)
1133 scaling_factor = 1/min(frob_matrix, gersh_matrix)
1134 matrix_sign = matrix*scaling_factor
1135 ALLOCATE (tmp1(sz, sz))
1136 ALLOCATE (tmp2(sz, sz))
1137
1138 converged = .false.
1139 DO i = 1, 100
1140 CALL dgemm('N', 'N', sz, sz, sz, -1.0_dp, matrix_sign, sz, matrix_sign, sz, 0.0_dp, tmp1, sz)
1141
1142 ! check convergence (frob norm of what should be the identity matrix minus identity matrix)
1143 frob_matrix_base = dlange('F', sz, sz, tmp1, sz, work)
1144 DO j = 1, sz
1145 tmp1(j, j) = tmp1(j, j) + 1.0_dp
1146 END DO
1147 frob_matrix = dlange('F', sz, sz, tmp1, sz, work)
1148
1149 IF (sign_order == 2) THEN
1150 prefactor = 0.5_dp
1151 ! update the above to 3*I-X*X
1152 DO j = 1, sz
1153 tmp1(j, j) = tmp1(j, j) + 2.0_dp
1154 END DO
1155 ELSE IF (sign_order == 3) THEN
1156 tmp2(:, :) = tmp1
1157 tmp1 = tmp1*4.0_dp/3.0_dp
1158 DO j = 1, sz
1159 tmp1(j, j) = tmp1(j, j) + 8.0_dp/3.0_dp
1160 END DO
1161 CALL dgemm('N', 'N', sz, sz, sz, 1.0_dp, tmp2, sz, tmp2, sz, 1.0_dp, tmp1, sz)
1162 prefactor = 3.0_dp/8.0_dp
1163 ELSE
1164 cpabort("requested order is not implemented.")
1165 END IF
1166
1167 CALL dgemm('N', 'N', sz, sz, sz, prefactor, matrix_sign, sz, tmp1, sz, 0.0_dp, tmp2, sz)
1168 matrix_sign = tmp2
1169
1170 ! frob_matrix/frob_matrix_base < SQRT(threshold)
1171 IF (frob_matrix*frob_matrix < (threshold*frob_matrix_base*frob_matrix_base)) THEN
1172 WRITE (unit_nr, '(T6,A,1X,I6,1X,A,1X,I3,E12.3)') &
1173 "Submatrix", matrix_id, "final NS sign iter", i, frob_matrix/frob_matrix_base
1174 CALL m_flush(unit_nr)
1175 converged = .true.
1176 EXIT
1177 END IF
1178 END DO
1179
1180 IF (.NOT. converged) THEN
1181 cpabort("dense_matrix_sign_Newton_Schulz did not converge within 100 iterations")
1182 END IF
1183
1184 DEALLOCATE (tmp1)
1185 DEALLOCATE (tmp2)
1186
1187 CALL timestop(handle)
1188
1189 END SUBROUTINE dense_matrix_sign_newton_schulz
1190
1191! **************************************************************************************************
1192!> \brief Perform eigendecomposition of a dense matrix
1193!> \param sm ...
1194!> \param N ...
1195!> \param eigvals ...
1196!> \param eigvecs ...
1197!> \par History
1198!> 2020.05 Extracted from dense_matrix_sign_direct [Michael Lass]
1199!> \author Michael Lass, Robert Schade
1200! **************************************************************************************************
1201 SUBROUTINE eigdecomp(sm, N, eigvals, eigvecs)
1202 INTEGER, INTENT(IN) :: n
1203 REAL(kind=dp), INTENT(IN) :: sm(n, n)
1204 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
1205 INTENT(OUT) :: eigvals
1206 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
1207 INTENT(OUT) :: eigvecs
1208
1209 INTEGER :: info, liwork, lwork
1210 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1211 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
1212 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp
1213
1214 ALLOCATE (eigvecs(n, n), tmp(n, n))
1215 ALLOCATE (eigvals(n))
1216
1217 ! symmetrize sm
1218 eigvecs(:, :) = 0.5*(sm + transpose(sm))
1219
1220 ! probe optimal sizes for WORK and IWORK
1221 lwork = -1
1222 liwork = -1
1223 ALLOCATE (work(1))
1224 ALLOCATE (iwork(1))
1225 CALL dsyevd('V', 'U', n, eigvecs, n, eigvals, work, lwork, iwork, liwork, info)
1226 lwork = int(work(1))
1227 liwork = int(iwork(1))
1228 DEALLOCATE (iwork)
1229 DEALLOCATE (work)
1230
1231 ! calculate eigenvalues and eigenvectors
1232 ALLOCATE (work(lwork))
1233 ALLOCATE (iwork(liwork))
1234 CALL dsyevd('V', 'U', n, eigvecs, n, eigvals, work, lwork, iwork, liwork, info)
1235 DEALLOCATE (iwork)
1236 DEALLOCATE (work)
1237 IF (info /= 0) cpabort("dsyevd did not succeed")
1238
1239 DEALLOCATE (tmp)
1240 END SUBROUTINE eigdecomp
1241
1242! **************************************************************************************************
1243!> \brief Calculate the sign matrix from eigenvalues and eigenvectors of a matrix
1244!> \param sm_sign ...
1245!> \param eigvals ...
1246!> \param eigvecs ...
1247!> \param N ...
1248!> \param mu_correction ...
1249!> \par History
1250!> 2020.05 Extracted from dense_matrix_sign_direct [Michael Lass]
1251!> \author Michael Lass, Robert Schade
1252! **************************************************************************************************
1253 SUBROUTINE sign_from_eigdecomp(sm_sign, eigvals, eigvecs, N, mu_correction)
1254 INTEGER :: n
1255 REAL(kind=dp), INTENT(IN) :: eigvecs(n, n), eigvals(n)
1256 REAL(kind=dp), INTENT(INOUT) :: sm_sign(n, n)
1257 REAL(kind=dp), INTENT(IN) :: mu_correction
1258
1259 INTEGER :: i
1260 REAL(kind=dp) :: modified_eigval, tmp(n, n)
1261
1262 sm_sign = 0
1263 DO i = 1, n
1264 modified_eigval = eigvals(i) - mu_correction
1265 IF (modified_eigval > 0) THEN
1266 sm_sign(i, i) = 1.0
1267 ELSE IF (modified_eigval < 0) THEN
1268 sm_sign(i, i) = -1.0
1269 ELSE
1270 sm_sign(i, i) = 0.0
1271 END IF
1272 END DO
1273
1274 ! Create matrix with eigenvalues in {-1,0,1} and eigenvectors of sm:
1275 ! sm_sign = eigvecs * sm_sign * eigvecs.T
1276 CALL dgemm('N', 'N', n, n, n, 1.0_dp, eigvecs, n, sm_sign, n, 0.0_dp, tmp, n)
1277 CALL dgemm('N', 'T', n, n, n, 1.0_dp, tmp, n, eigvecs, n, 0.0_dp, sm_sign, n)
1278 END SUBROUTINE sign_from_eigdecomp
1279
1280! **************************************************************************************************
1281!> \brief Compute partial trace of a matrix from its eigenvalues and eigenvectors
1282!> \param eigvals ...
1283!> \param eigvecs ...
1284!> \param firstcol ...
1285!> \param lastcol ...
1286!> \param mu_correction ...
1287!> \return ...
1288!> \par History
1289!> 2020.05 Created [Michael Lass]
1290!> \author Michael Lass
1291! **************************************************************************************************
1292 FUNCTION trace_from_eigdecomp(eigvals, eigvecs, firstcol, lastcol, mu_correction) RESULT(trace)
1293 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
1294 INTENT(IN) :: eigvals
1295 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
1296 INTENT(IN) :: eigvecs
1297 INTEGER, INTENT(IN) :: firstcol, lastcol
1298 REAL(kind=dp), INTENT(IN) :: mu_correction
1299 REAL(kind=dp) :: trace
1300
1301 INTEGER :: i, j, sm_size
1302 REAL(kind=dp) :: modified_eigval, tmpsum
1303 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: mapped_eigvals
1304
1305 sm_size = SIZE(eigvals)
1306 ALLOCATE (mapped_eigvals(sm_size))
1307
1308 DO i = 1, sm_size
1309 modified_eigval = eigvals(i) - mu_correction
1310 IF (modified_eigval > 0) THEN
1311 mapped_eigvals(i) = 1.0
1312 ELSE IF (modified_eigval < 0) THEN
1313 mapped_eigvals(i) = -1.0
1314 ELSE
1315 mapped_eigvals(i) = 0.0
1316 END IF
1317 END DO
1318
1319 trace = 0.0_dp
1320 DO i = firstcol, lastcol
1321 tmpsum = 0.0_dp
1322 DO j = 1, sm_size
1323 tmpsum = tmpsum + (eigvecs(i, j)*mapped_eigvals(j)*eigvecs(i, j))
1324 END DO
1325 trace = trace - 0.5_dp*tmpsum + 0.5_dp
1326 END DO
1327 END FUNCTION trace_from_eigdecomp
1328
1329! **************************************************************************************************
1330!> \brief Calculate the sign matrix by direct calculation of all eigenvalues and eigenvectors
1331!> \param sm_sign ...
1332!> \param sm ...
1333!> \param N ...
1334!> \par History
1335!> 2020.02 Created [Michael Lass, Robert Schade]
1336!> 2020.05 Extracted eigdecomp and sign_from_eigdecomp [Michael Lass]
1337!> \author Michael Lass, Robert Schade
1338! **************************************************************************************************
1339 SUBROUTINE dense_matrix_sign_direct(sm_sign, sm, N)
1340 INTEGER, INTENT(IN) :: n
1341 REAL(kind=dp), INTENT(IN) :: sm(n, n)
1342 REAL(kind=dp), INTENT(INOUT) :: sm_sign(n, n)
1343
1344 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvals
1345 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigvecs
1346
1347 CALL eigdecomp(sm, n, eigvals, eigvecs)
1348 CALL sign_from_eigdecomp(sm_sign, eigvals, eigvecs, n, 0.0_dp)
1349
1350 DEALLOCATE (eigvals, eigvecs)
1351 END SUBROUTINE dense_matrix_sign_direct
1352
1353! **************************************************************************************************
1354!> \brief Submatrix method
1355!> \param matrix_sign ...
1356!> \param matrix ...
1357!> \param threshold ...
1358!> \param sign_order ...
1359!> \param submatrix_sign_method ...
1360!> \par History
1361!> 2019.03 created [Robert Schade]
1362!> 2019.06 impl. submatrix method [Michael Lass]
1363!> \author Robert Schade, Michael Lass
1364! **************************************************************************************************
1365 SUBROUTINE matrix_sign_submatrix(matrix_sign, matrix, threshold, sign_order, submatrix_sign_method)
1366
1367 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_sign, matrix
1368 REAL(kind=dp), INTENT(IN) :: threshold
1369 INTEGER, INTENT(IN), OPTIONAL :: sign_order
1370 INTEGER, INTENT(IN) :: submatrix_sign_method
1371
1372 CHARACTER(LEN=*), PARAMETER :: routinen = 'matrix_sign_submatrix'
1373
1374 INTEGER :: handle, i, myrank, nblkcols, order, &
1375 sm_size, unit_nr
1376 INTEGER, ALLOCATABLE, DIMENSION(:) :: my_sms
1377 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: sm, sm_sign
1378 TYPE(cp_logger_type), POINTER :: logger
1379 TYPE(dbcsr_distribution_type) :: dist
1380 TYPE(submatrix_dissection_type) :: dissection
1381
1382 CALL timeset(routinen, handle)
1383
1384 ! print output on all ranks
1385 logger => cp_get_default_logger()
1386 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1387
1388 IF (PRESENT(sign_order)) THEN
1389 order = sign_order
1390 ELSE
1391 order = 2
1392 END IF
1393
1394 CALL dbcsr_get_info(matrix=matrix, nblkcols_total=nblkcols, distribution=dist)
1395 CALL dbcsr_distribution_get(dist=dist, mynode=myrank)
1396
1397 CALL dissection%init(matrix)
1398 CALL dissection%get_sm_ids_for_rank(myrank, my_sms)
1399
1400 !$OMP PARALLEL DEFAULT(OMP_DEFAULT_NONE_WITH_OOP) &
1401 !$OMP PRIVATE(sm, sm_sign, sm_size) &
1402 !$OMP SHARED(dissection, myrank, my_sms, order, submatrix_sign_method, threshold, unit_nr)
1403 !$OMP DO SCHEDULE(GUIDED)
1404 DO i = 1, SIZE(my_sms)
1405 WRITE (unit_nr, '(T3,A,1X,I4,1X,A,1X,I6)') "Rank", myrank, "processing submatrix", my_sms(i)
1406 CALL dissection%generate_submatrix(my_sms(i), sm)
1407 sm_size = SIZE(sm, 1)
1408 ALLOCATE (sm_sign(sm_size, sm_size))
1409 SELECT CASE (submatrix_sign_method)
1410 CASE (ls_scf_submatrix_sign_ns)
1411 CALL dense_matrix_sign_newton_schulz(sm_sign, sm, my_sms(i), threshold, order)
1412 CASE (ls_scf_submatrix_sign_direct, ls_scf_submatrix_sign_direct_muadj, ls_scf_submatrix_sign_direct_muadj_lowmem)
1413 CALL dense_matrix_sign_direct(sm_sign, sm, sm_size)
1414 CASE DEFAULT
1415 cpabort("Unkown submatrix sign method.")
1416 END SELECT
1417 CALL dissection%copy_resultcol(my_sms(i), sm_sign)
1418 DEALLOCATE (sm, sm_sign)
1419 END DO
1420 !$OMP END DO
1421 !$OMP END PARALLEL
1422
1423 CALL dissection%communicate_results(matrix_sign)
1424 CALL dissection%final
1425
1426 CALL timestop(handle)
1427
1428 END SUBROUTINE matrix_sign_submatrix
1429
1430! **************************************************************************************************
1431!> \brief Submatrix method with internal adjustment of chemical potential
1432!> \param matrix_sign ...
1433!> \param matrix ...
1434!> \param mu ...
1435!> \param nelectron ...
1436!> \param threshold ...
1437!> \param variant ...
1438!> \par History
1439!> 2020.05 Created [Michael Lass]
1440!> \author Robert Schade, Michael Lass
1441! **************************************************************************************************
1442 SUBROUTINE matrix_sign_submatrix_mu_adjust(matrix_sign, matrix, mu, nelectron, threshold, variant)
1443
1444 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_sign, matrix
1445 REAL(kind=dp), INTENT(INOUT) :: mu
1446 INTEGER, INTENT(IN) :: nelectron
1447 REAL(kind=dp), INTENT(IN) :: threshold
1448 INTEGER, INTENT(IN) :: variant
1449
1450 CHARACTER(LEN=*), PARAMETER :: routinen = 'matrix_sign_submatrix_mu_adjust'
1451 REAL(kind=dp), PARAMETER :: initial_increment = 0.01_dp
1452
1453 INTEGER :: handle, i, j, myrank, nblkcols, &
1454 sm_firstcol, sm_lastcol, sm_size, &
1455 unit_nr
1456 INTEGER, ALLOCATABLE, DIMENSION(:) :: my_sms
1457 LOGICAL :: has_mu_high, has_mu_low
1458 REAL(kind=dp) :: increment, mu_high, mu_low, new_mu, trace
1459 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: sm, sm_sign, tmp
1460 TYPE(cp_logger_type), POINTER :: logger
1461 TYPE(dbcsr_distribution_type) :: dist
1462 TYPE(eigbuf), ALLOCATABLE, DIMENSION(:) :: eigbufs
1463 TYPE(mp_comm_type) :: group
1464 TYPE(submatrix_dissection_type) :: dissection
1465
1466 CALL timeset(routinen, handle)
1467
1468 ! print output on all ranks
1469 logger => cp_get_default_logger()
1470 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1471
1472 CALL dbcsr_get_info(matrix=matrix, nblkcols_total=nblkcols, distribution=dist, group=group)
1473 CALL dbcsr_distribution_get(dist=dist, mynode=myrank)
1474
1475 CALL dissection%init(matrix)
1476 CALL dissection%get_sm_ids_for_rank(myrank, my_sms)
1477
1478 ALLOCATE (eigbufs(SIZE(my_sms)))
1479
1480 trace = 0.0_dp
1481
1482 !$OMP PARALLEL DEFAULT(OMP_DEFAULT_NONE_WITH_OOP) &
1483 !$OMP PRIVATE(sm, sm_sign, sm_size, sm_firstcol, sm_lastcol, j, tmp) &
1484 !$OMP SHARED(dissection, myrank, my_sms, unit_nr, eigbufs, threshold, variant) &
1485 !$OMP REDUCTION(+:trace)
1486 !$OMP DO SCHEDULE(GUIDED)
1487 DO i = 1, SIZE(my_sms)
1488 CALL dissection%generate_submatrix(my_sms(i), sm)
1489 sm_size = SIZE(sm, 1)
1490 WRITE (unit_nr, *) "Rank", myrank, "processing submatrix", my_sms(i), "size", sm_size
1491
1492 CALL dissection%get_relevant_sm_columns(my_sms(i), sm_firstcol, sm_lastcol)
1493
1494 IF (variant == ls_scf_submatrix_sign_direct_muadj) THEN
1495 ! Store all eigenvectors in buffer. We will use it to compute sm_sign at the end.
1496 CALL eigdecomp(sm, sm_size, eigvals=eigbufs(i)%eigvals, eigvecs=eigbufs(i)%eigvecs)
1497 ELSE
1498 ! Only store eigenvectors that are required for mu adjustment.
1499 ! Calculate sm_sign right away in the hope that mu is already correct.
1500 CALL eigdecomp(sm, sm_size, eigvals=eigbufs(i)%eigvals, eigvecs=tmp)
1501 ALLOCATE (eigbufs(i)%eigvecs(sm_firstcol:sm_lastcol, 1:sm_size))
1502 eigbufs(i)%eigvecs(:, :) = tmp(sm_firstcol:sm_lastcol, 1:sm_size)
1503
1504 ALLOCATE (sm_sign(sm_size, sm_size))
1505 CALL sign_from_eigdecomp(sm_sign, eigbufs(i)%eigvals, tmp, sm_size, 0.0_dp)
1506 CALL dissection%copy_resultcol(my_sms(i), sm_sign)
1507 DEALLOCATE (sm_sign, tmp)
1508 END IF
1509
1510 DEALLOCATE (sm)
1511 trace = trace + trace_from_eigdecomp(eigbufs(i)%eigvals, eigbufs(i)%eigvecs, sm_firstcol, sm_lastcol, 0.0_dp)
1512 END DO
1513 !$OMP END DO
1514 !$OMP END PARALLEL
1515
1516 has_mu_low = .false.
1517 has_mu_high = .false.
1518 increment = initial_increment
1519 new_mu = mu
1520 DO i = 1, 30
1521 CALL group%sum(trace)
1522 IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,1X,F13.9,1X,F15.9)') &
1523 "Density matrix: mu, trace error: ", new_mu, trace - nelectron
1524 IF (abs(trace - nelectron) < 0.5_dp) EXIT
1525 IF (trace < nelectron) THEN
1526 mu_low = new_mu
1527 new_mu = new_mu + increment
1528 has_mu_low = .true.
1529 increment = increment*2
1530 ELSE
1531 mu_high = new_mu
1532 new_mu = new_mu - increment
1533 has_mu_high = .true.
1534 increment = increment*2
1535 END IF
1536
1537 IF (has_mu_low .AND. has_mu_high) THEN
1538 new_mu = (mu_low + mu_high)/2
1539 IF (abs(mu_high - mu_low) < threshold) EXIT
1540 END IF
1541
1542 trace = 0
1543 !$OMP PARALLEL DEFAULT(OMP_DEFAULT_NONE_WITH_OOP) &
1544 !$OMP PRIVATE(i, sm_sign, tmp, sm_size, sm_firstcol, sm_lastcol) &
1545 !$OMP SHARED(dissection, my_sms, unit_nr, eigbufs, mu, new_mu, nelectron) &
1546 !$OMP REDUCTION(+:trace)
1547 !$OMP DO SCHEDULE(GUIDED)
1548 DO j = 1, SIZE(my_sms)
1549 sm_size = SIZE(eigbufs(j)%eigvals)
1550 CALL dissection%get_relevant_sm_columns(my_sms(j), sm_firstcol, sm_lastcol)
1551 trace = trace + trace_from_eigdecomp(eigbufs(j)%eigvals, eigbufs(j)%eigvecs, sm_firstcol, sm_lastcol, new_mu - mu)
1552 END DO
1553 !$OMP END DO
1554 !$OMP END PARALLEL
1555 END DO
1556
1557 ! Finalize sign matrix from eigendecompositions if we kept all eigenvectors
1558 IF (variant == ls_scf_submatrix_sign_direct_muadj) THEN
1559 !$OMP PARALLEL DEFAULT(OMP_DEFAULT_NONE_WITH_OOP) &
1560 !$OMP PRIVATE(sm, sm_sign, sm_size, sm_firstcol, sm_lastcol, j) &
1561 !$OMP SHARED(dissection, myrank, my_sms, unit_nr, eigbufs, mu, new_mu)
1562 !$OMP DO SCHEDULE(GUIDED)
1563 DO i = 1, SIZE(my_sms)
1564 WRITE (unit_nr, '(T3,A,1X,I4,1X,A,1X,I6)') "Rank", myrank, "finalizing submatrix", my_sms(i)
1565 sm_size = SIZE(eigbufs(i)%eigvals)
1566 ALLOCATE (sm_sign(sm_size, sm_size))
1567 CALL sign_from_eigdecomp(sm_sign, eigbufs(i)%eigvals, eigbufs(i)%eigvecs, sm_size, new_mu - mu)
1568 CALL dissection%copy_resultcol(my_sms(i), sm_sign)
1569 DEALLOCATE (sm_sign)
1570 END DO
1571 !$OMP END DO
1572 !$OMP END PARALLEL
1573 END IF
1574
1575 DEALLOCATE (eigbufs)
1576
1577 ! If we only stored parts of the eigenvectors and mu has changed, we need to recompute sm_sign
1578 IF ((variant == ls_scf_submatrix_sign_direct_muadj_lowmem) .AND. (mu /= new_mu)) THEN
1579 !$OMP PARALLEL DEFAULT(OMP_DEFAULT_NONE_WITH_OOP) &
1580 !$OMP PRIVATE(sm, sm_sign, sm_size, sm_firstcol, sm_lastcol, j) &
1581 !$OMP SHARED(dissection, myrank, my_sms, unit_nr, eigbufs, mu, new_mu)
1582 !$OMP DO SCHEDULE(GUIDED)
1583 DO i = 1, SIZE(my_sms)
1584 WRITE (unit_nr, '(T3,A,1X,I4,1X,A,1X,I6)') "Rank", myrank, "reprocessing submatrix", my_sms(i)
1585 CALL dissection%generate_submatrix(my_sms(i), sm)
1586 sm_size = SIZE(sm, 1)
1587 DO j = 1, sm_size
1588 sm(j, j) = sm(j, j) + mu - new_mu
1589 END DO
1590 ALLOCATE (sm_sign(sm_size, sm_size))
1591 CALL dense_matrix_sign_direct(sm_sign, sm, sm_size)
1592 CALL dissection%copy_resultcol(my_sms(i), sm_sign)
1593 DEALLOCATE (sm, sm_sign)
1594 END DO
1595 !$OMP END DO
1596 !$OMP END PARALLEL
1597 END IF
1598
1599 mu = new_mu
1600
1601 CALL dissection%communicate_results(matrix_sign)
1602 CALL dissection%final
1603
1604 CALL timestop(handle)
1605
1606 END SUBROUTINE matrix_sign_submatrix_mu_adjust
1607
1608! **************************************************************************************************
1609!> \brief compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations
1610!> the order of the algorithm should be 2..5, 3 or 5 is recommended
1611!> \param matrix_sqrt ...
1612!> \param matrix_sqrt_inv ...
1613!> \param matrix ...
1614!> \param threshold ...
1615!> \param order ...
1616!> \param eps_lanczos ...
1617!> \param max_iter_lanczos ...
1618!> \param symmetrize ...
1619!> \param converged ...
1620!> \param iounit ...
1621!> \par History
1622!> 2010.10 created [Joost VandeVondele]
1623!> \author Joost VandeVondele
1624! **************************************************************************************************
1625 SUBROUTINE matrix_sqrt_newton_schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, &
1626 eps_lanczos, max_iter_lanczos, symmetrize, converged, iounit)
1627 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_sqrt, matrix_sqrt_inv, matrix
1628 REAL(kind=dp), INTENT(IN) :: threshold
1629 INTEGER, INTENT(IN) :: order
1630 REAL(kind=dp), INTENT(IN) :: eps_lanczos
1631 INTEGER, INTENT(IN) :: max_iter_lanczos
1632 LOGICAL, OPTIONAL :: symmetrize, converged
1633 INTEGER, INTENT(IN), OPTIONAL :: iounit
1634
1635 CHARACTER(LEN=*), PARAMETER :: routinen = 'matrix_sqrt_Newton_Schulz'
1636
1637 INTEGER :: handle, i, unit_nr
1638 INTEGER(KIND=int_8) :: flop1, flop2, flop3, flop4, flop5
1639 LOGICAL :: arnoldi_converged, tsym
1640 REAL(kind=dp) :: a, b, c, conv, d, frob_matrix, &
1641 frob_matrix_base, gershgorin_norm, &
1642 max_ev, min_ev, oa, ob, oc, &
1643 occ_matrix, od, scaling, t1, t2
1644 TYPE(cp_logger_type), POINTER :: logger
1645 TYPE(dbcsr_type) :: tmp1, tmp2, tmp3
1646
1647 CALL timeset(routinen, handle)
1648
1649 IF (PRESENT(iounit)) THEN
1650 unit_nr = iounit
1651 ELSE
1652 logger => cp_get_default_logger()
1653 IF (logger%para_env%is_source()) THEN
1654 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1655 ELSE
1656 unit_nr = -1
1657 END IF
1658 END IF
1659
1660 IF (PRESENT(converged)) converged = .false.
1661 IF (PRESENT(symmetrize)) THEN
1662 tsym = symmetrize
1663 ELSE
1664 tsym = .true.
1665 END IF
1666
1667 ! for stability symmetry can not be assumed
1668 CALL dbcsr_create(tmp1, template=matrix, matrix_type=dbcsr_type_no_symmetry)
1669 CALL dbcsr_create(tmp2, template=matrix, matrix_type=dbcsr_type_no_symmetry)
1670 IF (order >= 4) THEN
1671 CALL dbcsr_create(tmp3, template=matrix, matrix_type=dbcsr_type_no_symmetry)
1672 END IF
1673
1674 CALL dbcsr_set(matrix_sqrt_inv, 0.0_dp)
1675 CALL dbcsr_add_on_diag(matrix_sqrt_inv, 1.0_dp)
1676 CALL dbcsr_filter(matrix_sqrt_inv, threshold)
1677 CALL dbcsr_copy(matrix_sqrt, matrix)
1678
1679 ! scale the matrix to get into the convergence range
1680 IF (order == 0) THEN
1681
1682 gershgorin_norm = dbcsr_gershgorin_norm(matrix_sqrt)
1683 frob_matrix = dbcsr_frobenius_norm(matrix_sqrt)
1684 scaling = 1.0_dp/min(frob_matrix, gershgorin_norm)
1685
1686 ELSE
1687
1688 ! scale the matrix to get into the convergence range
1689 CALL arnoldi_extremal(matrix_sqrt, max_ev, min_ev, threshold=eps_lanczos, &
1690 max_iter=max_iter_lanczos, converged=arnoldi_converged)
1691 IF (unit_nr > 0) THEN
1692 WRITE (unit_nr, *)
1693 WRITE (unit_nr, '(T6,A,1X,L1,A,E12.3)') "Lanczos converged: ", arnoldi_converged, " threshold:", eps_lanczos
1694 WRITE (unit_nr, '(T6,A,1X,E12.3,E12.3)') "Est. extremal eigenvalues:", max_ev, min_ev
1695 WRITE (unit_nr, '(T6,A,1X,E12.3)') "Est. condition number :", max_ev/max(min_ev, epsilon(min_ev))
1696 END IF
1697 ! conservatively assume we get a relatively large error (100*threshold_lanczos) in the estimates
1698 ! and adjust the scaling to be on the safe side
1699 scaling = 2.0_dp/(max_ev + min_ev + 100*eps_lanczos)
1700
1701 END IF
1702
1703 CALL dbcsr_scale(matrix_sqrt, scaling)
1704 CALL dbcsr_filter(matrix_sqrt, threshold)
1705 IF (unit_nr > 0) THEN
1706 WRITE (unit_nr, *)
1707 WRITE (unit_nr, *) "Order=", order
1708 END IF
1709
1710 DO i = 1, 100
1711
1712 t1 = m_walltime()
1713
1714 ! tmp1 = Zk * Yk - I
1715 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sqrt_inv, matrix_sqrt, 0.0_dp, tmp1, &
1716 filter_eps=threshold, flop=flop1)
1717 frob_matrix_base = dbcsr_frobenius_norm(tmp1)
1718 CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
1719
1720 ! check convergence (frob norm of what should be the identity matrix minus identity matrix)
1721 frob_matrix = dbcsr_frobenius_norm(tmp1)
1722
1723 flop4 = 0; flop5 = 0
1724 SELECT CASE (order)
1725 CASE (0, 2)
1726 ! update the above to 0.5*(3*I-Zk*Yk)
1727 CALL dbcsr_add_on_diag(tmp1, -2.0_dp)
1728 CALL dbcsr_scale(tmp1, -0.5_dp)
1729 CASE (3)
1730 ! tmp2 = tmp1 ** 2
1731 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp1, tmp1, 0.0_dp, tmp2, &
1732 filter_eps=threshold, flop=flop4)
1733 ! tmp1 = 1/16 * (16*I-8*tmp1+6*tmp1**2-5*tmp1**3)
1734 CALL dbcsr_add(tmp1, tmp2, -4.0_dp, 3.0_dp)
1735 CALL dbcsr_add_on_diag(tmp1, 8.0_dp)
1736 CALL dbcsr_scale(tmp1, 0.125_dp)
1737 CASE (4) ! as expensive as case(5), so little need to use it
1738 ! tmp2 = tmp1 ** 2
1739 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp1, tmp1, 0.0_dp, tmp2, &
1740 filter_eps=threshold, flop=flop4)
1741 ! tmp3 = tmp2 * tmp1
1742 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp2, tmp1, 0.0_dp, tmp3, &
1743 filter_eps=threshold, flop=flop5)
1744 CALL dbcsr_scale(tmp1, -8.0_dp)
1745 CALL dbcsr_add_on_diag(tmp1, 16.0_dp)
1746 CALL dbcsr_add(tmp1, tmp2, 1.0_dp, 6.0_dp)
1747 CALL dbcsr_add(tmp1, tmp3, 1.0_dp, -5.0_dp)
1748 CALL dbcsr_scale(tmp1, 1/16.0_dp)
1749 CASE (5)
1750 ! Knuth's reformulation to evaluate the polynomial of 4th degree in 2 multiplications
1751 ! p = y4+A*y3+B*y2+C*y+D
1752 ! z := y * (y+a); P := (z+y+b) * (z+c) + d.
1753 ! a=(A-1)/2 ; b=B*(a+1)-C-a*(a+1)*(a+1)
1754 ! c=B-b-a*(a+1)
1755 ! d=D-bc
1756 oa = -40.0_dp/35.0_dp
1757 ob = 48.0_dp/35.0_dp
1758 oc = -64.0_dp/35.0_dp
1759 od = 128.0_dp/35.0_dp
1760 a = (oa - 1)/2
1761 b = ob*(a + 1) - oc - a*(a + 1)**2
1762 c = ob - b - a*(a + 1)
1763 d = od - b*c
1764 ! tmp2 = tmp1 ** 2 + a * tmp1
1765 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp1, tmp1, 0.0_dp, tmp2, &
1766 filter_eps=threshold, flop=flop4)
1767 CALL dbcsr_add(tmp2, tmp1, 1.0_dp, a)
1768 ! tmp3 = tmp2 + tmp1 + b
1769 CALL dbcsr_copy(tmp3, tmp2)
1770 CALL dbcsr_add(tmp3, tmp1, 1.0_dp, 1.0_dp)
1771 CALL dbcsr_add_on_diag(tmp3, b)
1772 ! tmp2 = tmp2 + c
1773 CALL dbcsr_add_on_diag(tmp2, c)
1774 ! tmp1 = tmp2 * tmp3
1775 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp2, tmp3, 0.0_dp, tmp1, &
1776 filter_eps=threshold, flop=flop5)
1777 ! tmp1 = tmp1 + d
1778 CALL dbcsr_add_on_diag(tmp1, d)
1779 ! final scale
1780 CALL dbcsr_scale(tmp1, 35.0_dp/128.0_dp)
1781 CASE DEFAULT
1782 cpabort("Illegal order value")
1783 END SELECT
1784
1785 ! tmp2 = Yk * tmp1 = Y(k+1)
1786 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sqrt, tmp1, 0.0_dp, tmp2, &
1787 filter_eps=threshold, flop=flop2)
1788 ! CALL dbcsr_filter(tmp2,threshold)
1789 CALL dbcsr_copy(matrix_sqrt, tmp2)
1790
1791 ! tmp2 = tmp1 * Zk = Z(k+1)
1792 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp1, matrix_sqrt_inv, 0.0_dp, tmp2, &
1793 filter_eps=threshold, flop=flop3)
1794 ! CALL dbcsr_filter(tmp2,threshold)
1795 CALL dbcsr_copy(matrix_sqrt_inv, tmp2)
1796
1797 occ_matrix = dbcsr_get_occupation(matrix_sqrt_inv)
1798
1799 ! done iterating
1800 t2 = m_walltime()
1801
1802 conv = frob_matrix/frob_matrix_base
1803
1804 IF (unit_nr > 0) THEN
1805 WRITE (unit_nr, '(T6,A,1X,I3,1X,F10.8,E12.3,F12.3,F13.3)') "NS sqrt iter ", i, occ_matrix, &
1806 conv, t2 - t1, &
1807 (flop1 + flop2 + flop3 + flop4 + flop5)/(1.0e6_dp*max(0.001_dp, t2 - t1))
1808 CALL m_flush(unit_nr)
1809 END IF
1810
1811 IF (abnormal_value(conv)) THEN
1812 cpabort("conv is an abnormal value (NaN/Inf).")
1813 END IF
1814
1815 ! conv < SQRT(threshold)
1816 IF ((conv*conv) < threshold) THEN
1817 IF (PRESENT(converged)) converged = .true.
1818 EXIT
1819 END IF
1820
1821 END DO
1822
1823 ! symmetrize the matrices as this is not guaranteed by the algorithm
1824 IF (tsym) THEN
1825 IF (unit_nr > 0) THEN
1826 WRITE (unit_nr, '(T6,A20)') "Symmetrizing Results"
1827 END IF
1828 CALL dbcsr_transposed(tmp1, matrix_sqrt_inv)
1829 CALL dbcsr_add(matrix_sqrt_inv, tmp1, 0.5_dp, 0.5_dp)
1830 CALL dbcsr_transposed(tmp1, matrix_sqrt)
1831 CALL dbcsr_add(matrix_sqrt, tmp1, 0.5_dp, 0.5_dp)
1832 END IF
1833
1834 ! this check is not really needed
1835 CALL dbcsr_multiply("N", "N", +1.0_dp, matrix_sqrt_inv, matrix_sqrt, 0.0_dp, tmp1, &
1836 filter_eps=threshold)
1837 frob_matrix_base = dbcsr_frobenius_norm(tmp1)
1838 CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
1839 frob_matrix = dbcsr_frobenius_norm(tmp1)
1840 occ_matrix = dbcsr_get_occupation(matrix_sqrt_inv)
1841 IF (unit_nr > 0) THEN
1842 WRITE (unit_nr, '(T6,A,1X,I3,1X,F10.8,E12.3)') "Final NS sqrt iter ", i, occ_matrix, &
1843 frob_matrix/frob_matrix_base
1844 WRITE (unit_nr, '()')
1845 CALL m_flush(unit_nr)
1846 END IF
1847
1848 ! scale to proper end results
1849 CALL dbcsr_scale(matrix_sqrt, 1/sqrt(scaling))
1850 CALL dbcsr_scale(matrix_sqrt_inv, sqrt(scaling))
1851
1852 CALL dbcsr_release(tmp1)
1853 CALL dbcsr_release(tmp2)
1854 IF (order >= 4) THEN
1855 CALL dbcsr_release(tmp3)
1856 END IF
1857
1858 CALL timestop(handle)
1859
1860 END SUBROUTINE matrix_sqrt_newton_schulz
1861
1862! **************************************************************************************************
1863!> \brief compute the sqrt of a matrix via the general algorithm for the p-th root of Richters et al.
1864!> Commun. Comput. Phys., 25 (2019), pp. 564-585.
1865!> \param matrix_sqrt ...
1866!> \param matrix_sqrt_inv ...
1867!> \param matrix ...
1868!> \param threshold ...
1869!> \param order ...
1870!> \param eps_lanczos ...
1871!> \param max_iter_lanczos ...
1872!> \param symmetrize ...
1873!> \param converged ...
1874!> \par History
1875!> 2019.04 created [Robert Schade]
1876!> \author Robert Schade
1877! **************************************************************************************************
1878 SUBROUTINE matrix_sqrt_proot(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, &
1879 eps_lanczos, max_iter_lanczos, symmetrize, converged)
1880 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_sqrt, matrix_sqrt_inv, matrix
1881 REAL(kind=dp), INTENT(IN) :: threshold
1882 INTEGER, INTENT(IN) :: order
1883 REAL(kind=dp), INTENT(IN) :: eps_lanczos
1884 INTEGER, INTENT(IN) :: max_iter_lanczos
1885 LOGICAL, OPTIONAL :: symmetrize, converged
1886
1887 CHARACTER(LEN=*), PARAMETER :: routinen = 'matrix_sqrt_proot'
1888
1889 INTEGER :: choose, handle, i, ii, j, unit_nr
1890 INTEGER(KIND=int_8) :: f, flop1, flop2, flop3, flop4, flop5
1891 LOGICAL :: arnoldi_converged, test, tsym
1892 REAL(kind=dp) :: conv, frob_matrix, frob_matrix_base, &
1893 max_ev, min_ev, occ_matrix, scaling, &
1894 t1, t2
1895 TYPE(cp_logger_type), POINTER :: logger
1896 TYPE(dbcsr_type) :: bk2a, matrixs, rmat, tmp1, tmp2, tmp3
1897
1898 CALL cite_reference(richters2018)
1899
1900 test = .false.
1901
1902 CALL timeset(routinen, handle)
1903
1904 logger => cp_get_default_logger()
1905 IF (logger%para_env%is_source()) THEN
1906 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1907 ELSE
1908 unit_nr = -1
1909 END IF
1910
1911 IF (PRESENT(converged)) converged = .false.
1912 IF (PRESENT(symmetrize)) THEN
1913 tsym = symmetrize
1914 ELSE
1915 tsym = .true.
1916 END IF
1917
1918 ! for stability symmetry can not be assumed
1919 CALL dbcsr_create(tmp1, template=matrix, matrix_type=dbcsr_type_no_symmetry)
1920 CALL dbcsr_create(tmp2, template=matrix, matrix_type=dbcsr_type_no_symmetry)
1921 CALL dbcsr_create(tmp3, template=matrix, matrix_type=dbcsr_type_no_symmetry)
1922 CALL dbcsr_create(rmat, template=matrix, matrix_type=dbcsr_type_no_symmetry)
1923 CALL dbcsr_create(matrixs, template=matrix, matrix_type=dbcsr_type_no_symmetry)
1924
1925 CALL dbcsr_copy(matrixs, matrix)
1926 IF (1 == 1) THEN
1927 ! scale the matrix to get into the convergence range
1928 CALL arnoldi_extremal(matrixs, max_ev, min_ev, threshold=eps_lanczos, &
1929 max_iter=max_iter_lanczos, converged=arnoldi_converged)
1930 IF (unit_nr > 0) THEN
1931 WRITE (unit_nr, *)
1932 WRITE (unit_nr, '(T6,A,1X,L1,A,E12.3)') "Lanczos converged: ", arnoldi_converged, " threshold:", eps_lanczos
1933 WRITE (unit_nr, '(T6,A,1X,E12.3,E12.3)') "Est. extremal eigenvalues:", max_ev, min_ev
1934 WRITE (unit_nr, '(T6,A,1X,E12.3)') "Est. condition number :", max_ev/max(min_ev, epsilon(min_ev))
1935 END IF
1936 ! conservatively assume we get a relatively large error (100*threshold_lanczos) in the estimates
1937 ! and adjust the scaling to be on the safe side
1938 scaling = 2.0_dp/(max_ev + min_ev + 100*eps_lanczos)
1939 CALL dbcsr_scale(matrixs, scaling)
1940 CALL dbcsr_filter(matrixs, threshold)
1941 ELSE
1942 scaling = 1.0_dp
1943 END IF
1944
1945 CALL dbcsr_set(matrix_sqrt_inv, 0.0_dp)
1946 CALL dbcsr_add_on_diag(matrix_sqrt_inv, 1.0_dp)
1947 !CALL dbcsr_filter(matrix_sqrt_inv, threshold)
1948
1949 IF (unit_nr > 0) THEN
1950 WRITE (unit_nr, *)
1951 WRITE (unit_nr, *) "Order=", order
1952 END IF
1953
1954 DO i = 1, 100
1955
1956 t1 = m_walltime()
1957 IF (1 == 1) THEN
1958 !build R=1-A B_K^2
1959 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sqrt_inv, matrix_sqrt_inv, 0.0_dp, tmp1, &
1960 filter_eps=threshold, flop=flop1)
1961 CALL dbcsr_multiply("N", "N", 1.0_dp, matrixs, tmp1, 0.0_dp, rmat, &
1962 filter_eps=threshold, flop=flop2)
1963 CALL dbcsr_scale(rmat, -1.0_dp)
1964 CALL dbcsr_add_on_diag(rmat, 1.0_dp)
1965
1966 flop4 = 0; flop5 = 0
1967 CALL dbcsr_set(tmp1, 0.0_dp)
1968 CALL dbcsr_add_on_diag(tmp1, 2.0_dp)
1969
1970 flop3 = 0
1971
1972 DO j = 2, order
1973 IF (j == 2) THEN
1974 CALL dbcsr_copy(tmp2, rmat)
1975 ELSE
1976 f = 0
1977 CALL dbcsr_copy(tmp3, tmp2)
1978 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp3, rmat, 0.0_dp, tmp2, &
1979 filter_eps=threshold, flop=f)
1980 flop3 = flop3 + f
1981 END IF
1982 CALL dbcsr_add(tmp1, tmp2, 1.0_dp, 1.0_dp)
1983 END DO
1984 ELSE
1985 CALL dbcsr_create(bk2a, template=matrix, matrix_type=dbcsr_type_no_symmetry)
1986 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sqrt_inv, matrixs, 0.0_dp, tmp3, &
1987 filter_eps=threshold, flop=flop1)
1988 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sqrt_inv, tmp3, 0.0_dp, bk2a, &
1989 filter_eps=threshold, flop=flop2)
1990 CALL dbcsr_copy(rmat, bk2a)
1991 CALL dbcsr_add_on_diag(rmat, -1.0_dp)
1992
1993 CALL dbcsr_set(tmp1, 0.0_dp)
1994 CALL dbcsr_add_on_diag(tmp1, 1.0_dp)
1995
1996 CALL dbcsr_set(tmp2, 0.0_dp)
1997 CALL dbcsr_add_on_diag(tmp2, 1.0_dp)
1998
1999 flop3 = 0
2000 DO j = 1, order
2001 !choose=factorial(order)/(factorial(j)*factorial(order-j))
2002 choose = product([(ii, ii=1, order)])/(product([(ii, ii=1, j)])*product([(ii, ii=1, order - j)]))
2003 CALL dbcsr_add(tmp1, tmp2, 1.0_dp, -1.0_dp*(-1)**j*choose)
2004 IF (j < order) THEN
2005 f = 0
2006 CALL dbcsr_copy(tmp3, tmp2)
2007 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp3, bk2a, 0.0_dp, tmp2, &
2008 filter_eps=threshold, flop=f)
2009 flop3 = flop3 + f
2010 END IF
2011 END DO
2012 CALL dbcsr_release(bk2a)
2013 END IF
2014
2015 CALL dbcsr_copy(tmp3, matrix_sqrt_inv)
2016 CALL dbcsr_multiply("N", "N", 0.5_dp, tmp3, tmp1, 0.0_dp, matrix_sqrt_inv, &
2017 filter_eps=threshold, flop=flop4)
2018
2019 occ_matrix = dbcsr_get_occupation(matrix_sqrt_inv)
2020
2021 ! done iterating
2022 t2 = m_walltime()
2023
2024 conv = dbcsr_frobenius_norm(rmat)
2025
2026 IF (unit_nr > 0) THEN
2027 WRITE (unit_nr, '(T6,A,1X,I3,1X,F10.8,E12.3,F12.3,F13.3)') "PROOT sqrt iter ", i, occ_matrix, &
2028 conv, t2 - t1, &
2029 (flop1 + flop2 + flop3 + flop4 + flop5)/(1.0e6_dp*max(0.001_dp, t2 - t1))
2030 CALL m_flush(unit_nr)
2031 END IF
2032
2033 IF (abnormal_value(conv)) THEN
2034 cpabort("conv is an abnormal value (NaN/Inf).")
2035 END IF
2036
2037 ! conv < SQRT(threshold)
2038 IF ((conv*conv) < threshold) THEN
2039 IF (PRESENT(converged)) converged = .true.
2040 EXIT
2041 END IF
2042
2043 END DO
2044
2045 ! scale to proper end results
2046 CALL dbcsr_scale(matrix_sqrt_inv, sqrt(scaling))
2047 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sqrt_inv, matrix, 0.0_dp, matrix_sqrt, &
2048 filter_eps=threshold, flop=flop5)
2049
2050 ! symmetrize the matrices as this is not guaranteed by the algorithm
2051 IF (tsym) THEN
2052 IF (unit_nr > 0) THEN
2053 WRITE (unit_nr, '(A20)') "SYMMETRIZING RESULTS"
2054 END IF
2055 CALL dbcsr_transposed(tmp1, matrix_sqrt_inv)
2056 CALL dbcsr_add(matrix_sqrt_inv, tmp1, 0.5_dp, 0.5_dp)
2057 CALL dbcsr_transposed(tmp1, matrix_sqrt)
2058 CALL dbcsr_add(matrix_sqrt, tmp1, 0.5_dp, 0.5_dp)
2059 END IF
2060
2061 ! this check is not really needed
2062 IF (test) THEN
2063 CALL dbcsr_multiply("N", "N", +1.0_dp, matrix_sqrt_inv, matrix_sqrt, 0.0_dp, tmp1, &
2064 filter_eps=threshold)
2065 frob_matrix_base = dbcsr_frobenius_norm(tmp1)
2066 CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
2067 frob_matrix = dbcsr_frobenius_norm(tmp1)
2068 occ_matrix = dbcsr_get_occupation(matrix_sqrt_inv)
2069 IF (unit_nr > 0) THEN
2070 WRITE (unit_nr, '(T6,A,1X,I3,1X,F10.8,E12.3)') "Final PROOT S^{-1/2} S^{1/2}-Eins error ", i, occ_matrix, &
2071 frob_matrix/frob_matrix_base
2072 WRITE (unit_nr, '()')
2073 CALL m_flush(unit_nr)
2074 END IF
2075
2076 ! this check is not really needed
2077 CALL dbcsr_multiply("N", "N", +1.0_dp, matrix_sqrt_inv, matrix_sqrt_inv, 0.0_dp, tmp2, &
2078 filter_eps=threshold)
2079 CALL dbcsr_multiply("N", "N", +1.0_dp, tmp2, matrix, 0.0_dp, tmp1, &
2080 filter_eps=threshold)
2081 frob_matrix_base = dbcsr_frobenius_norm(tmp1)
2082 CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
2083 frob_matrix = dbcsr_frobenius_norm(tmp1)
2084 occ_matrix = dbcsr_get_occupation(matrix_sqrt_inv)
2085 IF (unit_nr > 0) THEN
2086 WRITE (unit_nr, '(T6,A,1X,I3,1X,F10.8,E12.3)') "Final PROOT S^{-1/2} S^{-1/2} S-Eins error ", i, occ_matrix, &
2087 frob_matrix/frob_matrix_base
2088 WRITE (unit_nr, '()')
2089 CALL m_flush(unit_nr)
2090 END IF
2091 END IF
2092
2093 CALL dbcsr_release(tmp1)
2094 CALL dbcsr_release(tmp2)
2095 CALL dbcsr_release(tmp3)
2096 CALL dbcsr_release(rmat)
2097 CALL dbcsr_release(matrixs)
2098
2099 CALL timestop(handle)
2100 END SUBROUTINE matrix_sqrt_proot
2101
2102! **************************************************************************************************
2103!> \brief ...
2104!> \param matrix_exp ...
2105!> \param matrix ...
2106!> \param omega ...
2107!> \param alpha ...
2108!> \param threshold ...
2109! **************************************************************************************************
2110 SUBROUTINE matrix_exponential(matrix_exp, matrix, omega, alpha, threshold)
2111 ! compute matrix_exp=omega*exp(alpha*matrix)
2112 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_exp, matrix
2113 REAL(kind=dp), INTENT(IN) :: omega, alpha, threshold
2114
2115 CHARACTER(LEN=*), PARAMETER :: routinen = 'matrix_exponential'
2116 REAL(dp), PARAMETER :: one = 1.0_dp, toll = 1.e-17_dp, &
2117 zero = 0.0_dp
2118
2119 INTEGER :: handle, i, k, unit_nr
2120 REAL(dp) :: factorial, norm_c, norm_d, norm_scalar
2121 TYPE(cp_logger_type), POINTER :: logger
2122 TYPE(dbcsr_type) :: b, b_square, c, d, d_product
2123
2124 CALL timeset(routinen, handle)
2125
2126 logger => cp_get_default_logger()
2127 IF (logger%para_env%is_source()) THEN
2128 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
2129 ELSE
2130 unit_nr = -1
2131 END IF
2132
2133 ! Calculate the norm of the matrix alpha*matrix, and scale it until it is less than 1.0
2134 norm_scalar = abs(alpha)*dbcsr_frobenius_norm(matrix)
2135
2136 ! k=scaling parameter
2137 k = 1
2138 DO
2139 IF ((norm_scalar/2.0_dp**k) <= one) EXIT
2140 k = k + 1
2141 END DO
2142
2143 ! copy and scale the input matrix in matrix C and in matrix D
2144 CALL dbcsr_create(c, template=matrix, matrix_type=dbcsr_type_no_symmetry)
2145 CALL dbcsr_copy(c, matrix)
2146 CALL dbcsr_scale(c, alpha_scalar=alpha/2.0_dp**k)
2147
2148 CALL dbcsr_create(d, template=matrix, matrix_type=dbcsr_type_no_symmetry)
2149 CALL dbcsr_copy(d, c)
2150
2151 ! write(*,*)
2152 ! write(*,*)
2153 ! CALL dbcsr_print(D, variable_name="D")
2154
2155 ! set the B matrix as B=Identity+D
2156 CALL dbcsr_create(b, template=matrix, matrix_type=dbcsr_type_no_symmetry)
2157 CALL dbcsr_copy(b, d)
2158 CALL dbcsr_add_on_diag(b, alpha=one)
2159
2160 ! CALL dbcsr_print(B, variable_name="B")
2161
2162 ! Calculate the norm of C and moltiply by toll to be used as a threshold
2163 norm_c = toll*dbcsr_frobenius_norm(matrix)
2164
2165 ! iteration for the truncated taylor series expansion
2166 CALL dbcsr_create(d_product, template=matrix, matrix_type=dbcsr_type_no_symmetry)
2167 i = 1
2168 DO
2169 i = i + 1
2170 ! compute D_product=D*C
2171 CALL dbcsr_multiply("N", "N", one, d, c, &
2172 zero, d_product, filter_eps=threshold)
2173
2174 ! copy D_product in D
2175 CALL dbcsr_copy(d, d_product)
2176
2177 ! calculate B=B+D_product/fat(i)
2178 factorial = ifac(i)
2179 CALL dbcsr_add(b, d_product, one, factorial)
2180
2181 ! check for convergence using the norm of D (copy of the matrix D_product) and C
2182 norm_d = factorial*dbcsr_frobenius_norm(d)
2183 IF (norm_d < norm_c) EXIT
2184 END DO
2185
2186 ! start the k iteration for the squaring of the matrix
2187 CALL dbcsr_create(b_square, template=matrix, matrix_type=dbcsr_type_no_symmetry)
2188 DO i = 1, k
2189 !compute B_square=B*B
2190 CALL dbcsr_multiply("N", "N", one, b, b, &
2191 zero, b_square, filter_eps=threshold)
2192 ! copy Bsquare in B to iterate
2193 CALL dbcsr_copy(b, b_square)
2194 END DO
2195
2196 ! copy B_square in matrix_exp and
2197 CALL dbcsr_copy(matrix_exp, b_square)
2198
2199 ! scale matrix_exp by omega, matrix_exp=omega*B_square
2200 CALL dbcsr_scale(matrix_exp, alpha_scalar=omega)
2201 ! write(*,*) alpha,omega
2202
2203 CALL dbcsr_release(b)
2204 CALL dbcsr_release(c)
2205 CALL dbcsr_release(d)
2206 CALL dbcsr_release(d_product)
2207 CALL dbcsr_release(b_square)
2208
2209 CALL timestop(handle)
2210
2211 END SUBROUTINE matrix_exponential
2212
2213! **************************************************************************************************
2214!> \brief McWeeny purification of a matrix in the orthonormal basis
2215!> \param matrix_p Matrix to purify (needs to be almost idempotent already)
2216!> \param threshold Threshold used as filter_eps and convergence criteria
2217!> \param max_steps Max number of iterations
2218!> \par History
2219!> 2013.01 created [Florian Schiffmann]
2220!> 2014.07 slightly refactored [Ole Schuett]
2221!> \author Florian Schiffmann
2222! **************************************************************************************************
2223 SUBROUTINE purify_mcweeny_orth(matrix_p, threshold, max_steps)
2224 TYPE(dbcsr_type), DIMENSION(:) :: matrix_p
2225 REAL(KIND=dp) :: threshold
2226 INTEGER :: max_steps
2227
2228 CHARACTER(LEN=*), PARAMETER :: routineN = 'purify_mcweeny_orth'
2229
2230 INTEGER :: handle, i, ispin, unit_nr
2231 REAL(KIND=dp) :: frob_norm, trace
2232 TYPE(cp_logger_type), POINTER :: logger
2233 TYPE(dbcsr_type) :: matrix_pp, matrix_tmp
2234
2235 CALL timeset(routinen, handle)
2236 logger => cp_get_default_logger()
2237 IF (logger%para_env%is_source()) THEN
2238 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
2239 ELSE
2240 unit_nr = -1
2241 END IF
2242
2243 CALL dbcsr_create(matrix_pp, template=matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
2244 CALL dbcsr_create(matrix_tmp, template=matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
2245 CALL dbcsr_trace(matrix_p(1), trace)
2246
2247 DO ispin = 1, SIZE(matrix_p)
2248 DO i = 1, max_steps
2249 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p(ispin), matrix_p(ispin), &
2250 0.0_dp, matrix_pp, filter_eps=threshold)
2251
2252 ! test convergence
2253 CALL dbcsr_copy(matrix_tmp, matrix_pp)
2254 CALL dbcsr_add(matrix_tmp, matrix_p(ispin), 1.0_dp, -1.0_dp)
2255 frob_norm = dbcsr_frobenius_norm(matrix_tmp) ! tmp = PP - P
2256 IF (unit_nr > 0) WRITE (unit_nr, '(t3,a,f16.8)') "McWeeny: Deviation of idempotency", frob_norm
2257 IF (unit_nr > 0) CALL m_flush(unit_nr)
2258
2259 ! construct new P
2260 CALL dbcsr_copy(matrix_tmp, matrix_pp)
2261 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_pp, matrix_p(ispin), &
2262 3.0_dp, matrix_tmp, filter_eps=threshold)
2263 CALL dbcsr_copy(matrix_p(ispin), matrix_tmp) ! tmp = 3PP - 2PPP
2264
2265 ! frob_norm < SQRT(trace*threshold)
2266 IF (frob_norm*frob_norm < trace*threshold) EXIT
2267 END DO
2268 END DO
2269
2270 CALL dbcsr_release(matrix_pp)
2271 CALL dbcsr_release(matrix_tmp)
2272 CALL timestop(handle)
2273 END SUBROUTINE purify_mcweeny_orth
2274
2275! **************************************************************************************************
2276!> \brief McWeeny purification of a matrix in the non-orthonormal basis
2277!> \param matrix_p Matrix to purify (needs to be almost idempotent already)
2278!> \param matrix_s Overlap-Matrix
2279!> \param threshold Threshold used as filter_eps and convergence criteria
2280!> \param max_steps Max number of iterations
2281!> \par History
2282!> 2013.01 created [Florian Schiffmann]
2283!> 2014.07 slightly refactored [Ole Schuett]
2284!> \author Florian Schiffmann
2285! **************************************************************************************************
2286 SUBROUTINE purify_mcweeny_nonorth(matrix_p, matrix_s, threshold, max_steps)
2287 TYPE(dbcsr_type), DIMENSION(:) :: matrix_p
2288 TYPE(dbcsr_type) :: matrix_s
2289 REAL(KIND=dp) :: threshold
2290 INTEGER :: max_steps
2291
2292 CHARACTER(LEN=*), PARAMETER :: routineN = 'purify_mcweeny_nonorth'
2293
2294 INTEGER :: handle, i, ispin, unit_nr
2295 REAL(KIND=dp) :: frob_norm, trace
2296 TYPE(cp_logger_type), POINTER :: logger
2297 TYPE(dbcsr_type) :: matrix_ps, matrix_psp, matrix_test
2298
2299 CALL timeset(routinen, handle)
2300
2301 logger => cp_get_default_logger()
2302 IF (logger%para_env%is_source()) THEN
2303 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
2304 ELSE
2305 unit_nr = -1
2306 END IF
2307
2308 CALL dbcsr_create(matrix_ps, template=matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
2309 CALL dbcsr_create(matrix_psp, template=matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
2310 CALL dbcsr_create(matrix_test, template=matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
2311
2312 DO ispin = 1, SIZE(matrix_p)
2313 DO i = 1, max_steps
2314 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p(ispin), matrix_s, &
2315 0.0_dp, matrix_ps, filter_eps=threshold)
2316 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ps, matrix_p(ispin), &
2317 0.0_dp, matrix_psp, filter_eps=threshold)
2318 IF (i == 1) CALL dbcsr_trace(matrix_ps, trace)
2319
2320 ! test convergence
2321 CALL dbcsr_copy(matrix_test, matrix_psp)
2322 CALL dbcsr_add(matrix_test, matrix_p(ispin), 1.0_dp, -1.0_dp)
2323 frob_norm = dbcsr_frobenius_norm(matrix_test) ! test = PSP - P
2324 IF (unit_nr > 0) WRITE (unit_nr, '(t3,a,2f16.8)') "McWeeny: Deviation of idempotency", frob_norm
2325 IF (unit_nr > 0) CALL m_flush(unit_nr)
2326
2327 ! construct new P
2328 CALL dbcsr_copy(matrix_p(ispin), matrix_psp)
2329 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_ps, matrix_psp, &
2330 3.0_dp, matrix_p(ispin), filter_eps=threshold)
2331
2332 ! frob_norm < SQRT(trace*threshold)
2333 IF (frob_norm*frob_norm < trace*threshold) EXIT
2334 END DO
2335 END DO
2336
2337 CALL dbcsr_release(matrix_ps)
2338 CALL dbcsr_release(matrix_psp)
2339 CALL dbcsr_release(matrix_test)
2340 CALL timestop(handle)
2341 END SUBROUTINE purify_mcweeny_nonorth
2342
2343END MODULE iterate_matrix
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
iterate_matrix::purify_mcweeny
Definition iterate_matrix.F:58
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:16
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
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::richters2018
integer, save, public richters2018
Definition bibliography.F:36
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_transposed
subroutine, public dbcsr_transposed(transposed, normal, shallow_data_copy, transpose_distribution, use_distribution)
...
Definition cp_dbcsr_api.F:1242
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1188
cp_dbcsr_api::dbcsr_get_matrix_type
character function, public dbcsr_get_matrix_type(matrix)
...
Definition cp_dbcsr_api.F:861
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:525
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:658
cp_dbcsr_api::dbcsr_get_occupation
real(kind=dp) function, public dbcsr_get_occupation(matrix)
...
Definition cp_dbcsr_api.F:893
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_api::dbcsr_distribution_get
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Definition cp_dbcsr_api.F:572
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_set_diag
subroutine, public dbcsr_set_diag(matrix, diag)
Copies the diagonal elements from the given array into the given matrix.
Definition cp_dbcsr_contrib.F:585
cp_dbcsr_contrib::dbcsr_gershgorin_norm
real(dp) function, public dbcsr_gershgorin_norm(matrix)
Compute the gershgorin norm of a dbcsr matrix.
Definition cp_dbcsr_contrib.F:167
cp_dbcsr_contrib::dbcsr_get_diag
subroutine, public dbcsr_get_diag(matrix, diag)
Copies the diagonal elements from the given matrix into the given array.
Definition cp_dbcsr_contrib.F:544
cp_dbcsr_contrib::dbcsr_add_on_diag
subroutine, public dbcsr_add_on_diag(matrix, alpha)
Adds the given scalar to the diagonal of the matrix. Reserves any missing diagonal blocks.
Definition cp_dbcsr_contrib.F:331
cp_dbcsr_contrib::dbcsr_maxabs
real(dp) function, public dbcsr_maxabs(matrix)
Compute the maxabs norm of a dbcsr matrix.
Definition cp_dbcsr_contrib.F:94
cp_dbcsr_contrib::dbcsr_trace
subroutine, public dbcsr_trace(matrix, trace)
Computes the trace of the given matrix, also known as the sum of its diagonal elements.
Definition cp_dbcsr_contrib.F:412
cp_dbcsr_contrib::dbcsr_frobenius_norm
real(dp) function, public dbcsr_frobenius_norm(matrix)
Compute the frobenius norm of a dbcsr matrix.
Definition cp_dbcsr_contrib.F:127
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:576
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::ls_scf_submatrix_sign_direct_muadj
integer, parameter, public ls_scf_submatrix_sign_direct_muadj
Definition input_constants.F:1022
input_constants::ls_scf_submatrix_sign_ns
integer, parameter, public ls_scf_submatrix_sign_ns
Definition input_constants.F:1022
input_constants::ls_scf_submatrix_sign_direct_muadj_lowmem
integer, parameter, public ls_scf_submatrix_sign_direct_muadj_lowmem
Definition input_constants.F:1022
input_constants::ls_scf_submatrix_sign_direct
integer, parameter, public ls_scf_submatrix_sign_direct
Definition input_constants.F:1022
iterate_matrix
Routines useful for iterative matrix calculations.
Definition iterate_matrix.F:13
iterate_matrix::determinant
recursive subroutine, public determinant(matrix, det, threshold)
Computes the determinant of a symmetric positive definite matrix using the trace of the matrix logari...
Definition iterate_matrix.F:85
iterate_matrix::matrix_sign_newton_schulz
subroutine, public matrix_sign_newton_schulz(matrix_sign, matrix, threshold, sign_order, iounit)
compute the sign a matrix using Newton-Schulz iterations
Definition iterate_matrix.F:682
iterate_matrix::matrix_sign_submatrix
subroutine, public matrix_sign_submatrix(matrix_sign, matrix, threshold, sign_order, submatrix_sign_method)
Submatrix method.
Definition iterate_matrix.F:1366
iterate_matrix::invert_taylor
subroutine, public invert_taylor(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite diagonally dominant matrix
Definition iterate_matrix.F:287
iterate_matrix::matrix_exponential
subroutine, public matrix_exponential(matrix_exp, matrix, omega, alpha, threshold)
...
Definition iterate_matrix.F:2111
iterate_matrix::matrix_sign_submatrix_mu_adjust
subroutine, public matrix_sign_submatrix_mu_adjust(matrix_sign, matrix, mu, nelectron, threshold, variant)
Submatrix method with internal adjustment of chemical potential.
Definition iterate_matrix.F:1443
iterate_matrix::invert_hotelling
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
Definition iterate_matrix.F:473
iterate_matrix::matrix_sign_proot
subroutine, public matrix_sign_proot(matrix_sign, matrix, threshold, sign_order)
compute the sign a matrix using the general algorithm for the p-th root of Richters et al....
Definition iterate_matrix.F:1014
iterate_matrix::matrix_sqrt_newton_schulz
subroutine, public matrix_sqrt_newton_schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged, iounit)
compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations the...
Definition iterate_matrix.F:1627
iterate_matrix::matrix_sqrt_proot
subroutine, public matrix_sqrt_proot(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged)
compute the sqrt of a matrix via the general algorithm for the p-th root of Richters et al....
Definition iterate_matrix.F:1880
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::ifac
real(kind=dp), dimension(0:maxfac), parameter, public ifac
Definition mathconstants.F:49
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::abnormal_value
logical function, public abnormal_value(a)
determines if a value is not normal (e.g. for Inf and Nan) based on IO to work also under optimizatio...
Definition mathlib.F:159
matrix_exp
Routines for calculating a complex matrix exponential.
Definition matrix_exp.F:13
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
submatrix_dissection
Definition submatrix_dissection.F:40
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_comm_type
Definition message_passing.F:161
submatrix_dissection::submatrix_dissection_type
Definition submatrix_dissection.F:66