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