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