(git:ed6f26b)
Loading...
Searching...
No Matches
matrix_exp.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! **************************************************************************************************
9!> \brief Routines for calculating a complex matrix exponential.
10!> \author Florian Schiffmann (02.09)
11! **************************************************************************************************
12
13MODULE matrix_exp
14
15 USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add,&
16 cp_cfm_solve
17 USE cp_cfm_types, ONLY: cp_cfm_create,&
18 cp_cfm_release,&
19 cp_cfm_set_all,&
20 cp_cfm_to_cfm,&
21 cp_cfm_type
22 USE cp_fm_basic_linalg, ONLY: cp_complex_fm_gemm,&
23 cp_fm_scale_and_add,&
24 cp_fm_solve
25 USE cp_fm_struct, ONLY: cp_fm_struct_double,&
26 cp_fm_struct_release,&
27 cp_fm_struct_type
28 USE cp_fm_types, ONLY: cp_fm_create,&
29 cp_fm_get_info,&
30 cp_fm_release,&
31 cp_fm_set_all,&
32 cp_fm_to_fm,&
33 cp_fm_type
34 USE cp_log_handling, ONLY: cp_to_string
35 USE kinds, ONLY: dp
36 USE mathconstants, ONLY: fac,&
37 z_one,&
38 z_zero
39 USE message_passing, ONLY: mp_comm_type,&
40 mp_para_env_type
41 USE parallel_gemm_api, ONLY: parallel_gemm
42#include "./base/base_uses.f90"
43
44 IMPLICIT NONE
45
46 PRIVATE
47
48 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'matrix_exp'
49
50 PUBLIC :: taylor_only_imaginary, &
51 taylor_full_complex, &
52 exp_pade_full_complex, &
53 exp_pade_only_imaginary, &
54 get_nsquare_norder, &
55 arnoldi, exp_pade_real
56
57CONTAINS
58
59! **************************************************************************************************
60!> \brief specialized subroutine for purely imaginary matrix exponentials
61!> \param exp_H ...
62!> \param im_matrix ...
63!> \param nsquare ...
64!> \param ntaylor ...
65!> \author Florian Schiffmann (02.09)
66! **************************************************************************************************
67
68 SUBROUTINE taylor_only_imaginary(exp_H, im_matrix, nsquare, ntaylor)
69 TYPE(cp_fm_type), DIMENSION(2) :: exp_h
70 TYPE(cp_fm_type), INTENT(IN) :: im_matrix
71 INTEGER, INTENT(in) :: nsquare, ntaylor
72
73 CHARACTER(len=*), PARAMETER :: routinen = 'taylor_only_imaginary'
74 REAL(kind=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
75
76 INTEGER :: handle, i, ndim, nloop
77 REAL(kind=dp) :: square_fac, tfac, tmp
78 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
79 POINTER :: local_data_im
80 TYPE(cp_fm_type) :: t1, t2, tres_im, tres_re
81
82 CALL timeset(routinen, handle)
83
84 CALL cp_fm_get_info(im_matrix, local_data=local_data_im)
85 ndim = im_matrix%matrix_struct%nrow_global
86
87 square_fac = 1.0_dp/(2.0_dp**real(nsquare, dp))
88! CALL cp_fm_scale(square_fac,im_matrix)
89 CALL cp_fm_create(t1, &
90 matrix_struct=im_matrix%matrix_struct, &
91 name="T1")
92
93 CALL cp_fm_create(t2, &
94 matrix_struct=t1%matrix_struct, &
95 name="T2")
96 CALL cp_fm_create(tres_im, &
97 matrix_struct=t1%matrix_struct, &
98 name="T3")
99 CALL cp_fm_create(tres_re, &
100 matrix_struct=t1%matrix_struct, &
101 name="Tres")
102 tmp = 1.0_dp
103
104 CALL cp_fm_set_all(tres_re, zero, one)
105 CALL cp_fm_set_all(tres_im, zero, zero)
106 CALL cp_fm_set_all(t1, zero, one)
107
108 tfac = one
109 nloop = ceiling(real(ntaylor, dp)/2.0_dp)
110
111 DO i = 1, nloop
112 tmp = tmp*(real(i, dp)*2.0_dp - 1.0_dp)
113 CALL parallel_gemm("N", "N", ndim, ndim, ndim, square_fac, im_matrix, t1, zero, &
114 ! CALL parallel_gemm("N","N",ndim,ndim,ndim,one,im_matrix,T1,zero,&
115 t2)
116 tfac = 1._dp/tmp
117 IF (mod(i, 2) == 0) tfac = -tfac
118 CALL cp_fm_scale_and_add(one, tres_im, tfac, t2)
119 tmp = tmp*real(i, dp)*2.0_dp
120 CALL parallel_gemm("N", "N", ndim, ndim, ndim, square_fac, im_matrix, t2, zero, &
121 ! CALL parallel_gemm("N","N",ndim,ndim,ndim,one,im_matrix,T2,zero,&
122 t1)
123 tfac = 1._dp/tmp
124 IF (mod(i, 2) == 1) tfac = -tfac
125 CALL cp_fm_scale_and_add(one, tres_re, tfac, t1)
126
127 END DO
128
129 IF (nsquare .GT. 0) THEN
130 DO i = 1, nsquare
131 CALL cp_complex_fm_gemm("N", "N", ndim, ndim, ndim, one, tres_re, tres_im, &
132 tres_re, tres_im, zero, exp_h(1), exp_h(2))
133
134 CALL cp_fm_to_fm(exp_h(1), tres_re)
135 CALL cp_fm_to_fm(exp_h(2), tres_im)
136 END DO
137 ELSE
138 CALL cp_fm_to_fm(tres_re, exp_h(1))
139 CALL cp_fm_to_fm(tres_im, exp_h(2))
140 END IF
141
142 CALL cp_fm_release(t1)
143 CALL cp_fm_release(t2)
144 CALL cp_fm_release(tres_re)
145 CALL cp_fm_release(tres_im)
146
147 CALL timestop(handle)
148
149 END SUBROUTINE taylor_only_imaginary
150
151! **************************************************************************************************
152!> \brief subroutine for general complex matrix exponentials
153!> on input a separate cp_fm_type for real and complex part
154!> on output a size 2 cp_fm_type, first element is the real part of
155!> the exponential second the imaginary
156!> \param exp_H ...
157!> \param re_part ...
158!> \param im_part ...
159!> \param nsquare ...
160!> \param ntaylor ...
161!> \author Florian Schiffmann (02.09)
162! **************************************************************************************************
163
164 SUBROUTINE taylor_full_complex(exp_H, re_part, im_part, nsquare, ntaylor)
165 TYPE(cp_fm_type), DIMENSION(2) :: exp_h
166 TYPE(cp_fm_type), INTENT(IN) :: re_part, im_part
167 INTEGER, INTENT(in) :: nsquare, ntaylor
168
169 CHARACTER(len=*), PARAMETER :: routinen = 'taylor_full_complex'
170
171 COMPLEX(KIND=dp) :: tfac
172 INTEGER :: handle, i, ndim
173 REAL(kind=dp) :: square_fac, tmp
174 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
175 POINTER :: local_data_im, local_data_re
176 TYPE(cp_cfm_type) :: t1, t2, t3, tres
177
178 CALL timeset(routinen, handle)
179 CALL cp_fm_get_info(re_part, local_data=local_data_re)
180 CALL cp_fm_get_info(im_part, local_data=local_data_im)
181 ndim = re_part%matrix_struct%nrow_global
182
183 CALL cp_cfm_create(t1, &
184 matrix_struct=re_part%matrix_struct, &
185 name="T1")
186
187 square_fac = 2.0_dp**real(nsquare, dp)
188
189 t1%local_data = cmplx(local_data_re/square_fac, local_data_im/square_fac, kind=dp)
190
191 CALL cp_cfm_create(t2, &
192 matrix_struct=t1%matrix_struct, &
193 name="T2")
194 CALL cp_cfm_create(t3, &
195 matrix_struct=t1%matrix_struct, &
196 name="T3")
197 CALL cp_cfm_create(tres, &
198 matrix_struct=t1%matrix_struct, &
199 name="Tres")
200 tmp = 1.0_dp
201 CALL cp_cfm_set_all(tres, z_zero, z_one)
202 CALL cp_cfm_set_all(t2, z_zero, z_one)
203 tfac = z_one
204
205 DO i = 1, ntaylor
206 tmp = tmp*real(i, dp)
207 CALL parallel_gemm("N", "N", ndim, ndim, ndim, z_one, t1, t2, z_zero, &
208 t3)
209 tfac = cmplx(1._dp/tmp, 0.0_dp, kind=dp)
210 CALL cp_cfm_scale_and_add(z_one, tres, tfac, t3)
211 CALL cp_cfm_to_cfm(t3, t2)
212 END DO
213
214 IF (nsquare .GT. 0) THEN
215 DO i = 1, nsquare
216 CALL parallel_gemm("N", "N", ndim, ndim, ndim, z_one, tres, tres, z_zero, &
217 t2)
218 CALL cp_cfm_to_cfm(t2, tres)
219 END DO
220 END IF
221
222 exp_h(1)%local_data = real(tres%local_data, kind=dp)
223 exp_h(2)%local_data = aimag(tres%local_data)
224
225 CALL cp_cfm_release(t1)
226 CALL cp_cfm_release(t2)
227 CALL cp_cfm_release(t3)
228 CALL cp_cfm_release(tres)
229 CALL timestop(handle)
230
231 END SUBROUTINE taylor_full_complex
232
233! **************************************************************************************************
234!> \brief optimization function for pade/taylor order and number of squaring steps
235!> \param norm ...
236!> \param nsquare ...
237!> \param norder ...
238!> \param eps_exp ...
239!> \param method ...
240!> \param do_emd ...
241!> \author Florian Schiffmann (02.09)
242! **************************************************************************************************
243 SUBROUTINE get_nsquare_norder(norm, nsquare, norder, eps_exp, method, do_emd)
244
245 REAL(dp), INTENT(in) :: norm
246 INTEGER, INTENT(out) :: nsquare, norder
247 REAL(dp), INTENT(in) :: eps_exp
248 INTEGER, INTENT(in) :: method
249 LOGICAL, INTENT(in) :: do_emd
250
251 INTEGER :: cost, i, iscale, orders(3), p, &
252 prev_cost, q
253 LOGICAL :: new_scale
254 REAL(dp) :: d, eval, myval, n, scaled, scalen
255
256 orders(:) = (/12, 12, 12/)
257 IF (method == 2) THEN
258 DO iscale = 0, 12
259 new_scale = .false.
260 eval = norm/(2.0_dp**real(iscale, dp))
261 DO q = 1, 12
262 DO p = max(1, q - 1), q
263 IF (p > q) EXIT
264 d = 1.0_dp
265 n = 1.0_dp
266 DO i = 1, q
267 IF (i .LE. p) scalen = fac(p + q - i)*fac(p)/(fac(p + q)*fac(i)*fac(p - i))
268 scaled = (-1.0)**i*fac(p + q - i)*fac(q)/(fac(p + q)*fac(i)*fac(q - i))
269 IF (i .LE. p) n = n + scalen*eval**i
270 d = d + scaled*eval**i
271 END DO
272 IF (abs((exp(norm) - (n/d)**(2.0_dp**iscale))/max(1.0_dp, exp(norm))) .LE. eps_exp) THEN
273 IF (do_emd) THEN
274 cost = iscale + q
275 prev_cost = orders(1) + orders(2)
276 ELSE
277 cost = iscale + ceiling(real(q, dp)/3.0_dp)
278 prev_cost = orders(1) + ceiling(real(orders(2), dp)/3.0_dp)
279 END IF
280 IF (cost .LT. prev_cost) THEN
281 orders(:) = (/iscale, q, p/)
282 myval = (n/d)**(2.0_dp**iscale)
283 END IF
284 new_scale = .true.
285 EXIT
286 END IF
287 END DO
288 IF (new_scale) EXIT
289 END DO
290 IF (iscale .GE. orders(1) + orders(2) .AND. new_scale) EXIT
291 END DO
292 ELSE IF (method == 1) THEN
293 q = 0
294 eval = norm
295 DO iscale = 0, 6
296 new_scale = .false.
297 IF (iscale .GE. 1) eval = norm/(2.0_dp**real(iscale, dp))
298 DO p = 1, 20
299 d = 1.0_dp
300 n = 1.0_dp
301 DO i = 1, p
302 scalen = 1.0_dp/fac(i)
303 n = n + scalen*(eval**real(i, dp))
304 END DO
305 IF (abs((exp(norm) - n**(2.0_dp**real(iscale, dp)))/max(1.0_dp, exp(norm))) .LE. eps_exp) THEN
306 IF (do_emd) THEN
307 cost = iscale + p
308 prev_cost = orders(1) + orders(2)
309 ELSE
310 cost = iscale + ceiling(real(p, dp)/3.0_dp)
311 prev_cost = orders(1) + ceiling(real(orders(2), dp)/3.0_dp)
312 END IF
313 IF (cost .LT. prev_cost) THEN
314 orders(:) = (/iscale, p, 0/)
315 myval = (n)**(2.0_dp**iscale)
316 END IF
317 new_scale = .true.
318 EXIT
319 END IF
320 END DO
321 IF (iscale .GE. orders(1) + orders(2) .AND. new_scale) EXIT
322 END DO
323 END IF
324
325 nsquare = orders(1)
326 norder = orders(2)
327
328 END SUBROUTINE get_nsquare_norder
329
330! **************************************************************************************************
331!> \brief exponential of a complex matrix,
332!> calculated using pade approximation together with scaling and squaring
333!> \param exp_H ...
334!> \param re_part ...
335!> \param im_part ...
336!> \param nsquare ...
337!> \param npade ...
338!> \author Florian Schiffmann (02.09)
339! **************************************************************************************************
340
341 SUBROUTINE exp_pade_full_complex(exp_H, re_part, im_part, nsquare, npade)
342 TYPE(cp_fm_type), DIMENSION(2) :: exp_h
343 TYPE(cp_fm_type), INTENT(IN) :: re_part, im_part
344 INTEGER, INTENT(in) :: nsquare, npade
345
346 CHARACTER(len=*), PARAMETER :: routinen = 'exp_pade_full_complex'
347
348 COMPLEX(KIND=dp) :: scaled, scalen
349 INTEGER :: handle, i, ldim, ndim, p, q
350 REAL(kind=dp) :: square_fac, tmp
351 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
352 POINTER :: local_data_im, local_data_re
353 TYPE(cp_cfm_type) :: dpq, fin_p, t1
354 TYPE(cp_cfm_type), DIMENSION(2) :: mult_p
355 TYPE(cp_cfm_type), TARGET :: npq, t2, tres
356
357 p = npade
358 q = npade
359
360 CALL timeset(routinen, handle)
361 CALL cp_fm_get_info(re_part, local_data=local_data_re, ncol_local=ldim, &
362 nrow_global=ndim)
363 CALL cp_fm_get_info(im_part, local_data=local_data_im)
364
365 CALL cp_cfm_create(dpq, &
366 matrix_struct=re_part%matrix_struct, &
367 name="Dpq")
368
369 square_fac = 2.0_dp**real(nsquare, dp)
370
371 CALL cp_cfm_create(t1, &
372 matrix_struct=dpq%matrix_struct, &
373 name="T1")
374
375 CALL cp_cfm_create(t2, &
376 matrix_struct=t1%matrix_struct, &
377 name="T2")
378 CALL cp_cfm_create(npq, &
379 matrix_struct=t1%matrix_struct, &
380 name="Npq")
381 CALL cp_cfm_create(tres, &
382 matrix_struct=t1%matrix_struct, &
383 name="Tres")
384
385 DO i = 1, ldim
386 t2%local_data(:, i) = cmplx(local_data_re(:, i)/square_fac, local_data_im(:, i)/square_fac, kind=dp)
387 END DO
388 CALL cp_cfm_to_cfm(t2, t1)
389 mult_p(1) = t2
390 mult_p(2) = tres
391 tmp = 1.0_dp
392 CALL cp_cfm_set_all(npq, z_zero, z_one)
393 CALL cp_cfm_set_all(dpq, z_zero, z_one)
394
395 CALL cp_cfm_scale_and_add(z_one, npq, z_one*0.5_dp, t2)
396 CALL cp_cfm_scale_and_add(z_one, dpq, -z_one*0.5_dp, t2)
397
398 IF (npade .GT. 2) THEN
399 DO i = 2, npade
400 IF (i .LE. p) scalen = cmplx(fac(p + q - i)*fac(p)/(fac(p + q)*fac(i)*fac(p - i)), 0.0_dp, kind=dp)
401 scaled = cmplx((-1.0_dp)**i*fac(p + q - i)*fac(q)/(fac(p + q)*fac(i)*fac(q - i)), 0.0_dp, kind=dp)
402 CALL parallel_gemm("N", "N", ndim, ndim, ndim, z_one, t1, mult_p(mod(i, 2) + 1), z_zero, &
403 mult_p(mod(i + 1, 2) + 1))
404 IF (i .LE. p) CALL cp_cfm_scale_and_add(z_one, npq, scalen, mult_p(mod(i + 1, 2) + 1))
405 IF (i .LE. q) CALL cp_cfm_scale_and_add(z_one, dpq, scaled, mult_p(mod(i + 1, 2) + 1))
406 END DO
407 END IF
408
409 CALL cp_cfm_solve(dpq, npq)
410
411 mult_p(2) = npq
412 mult_p(1) = tres
413 IF (nsquare .GT. 0) THEN
414 DO i = 1, nsquare
415 CALL parallel_gemm("N", "N", ndim, ndim, ndim, z_one, mult_p(mod(i, 2) + 1), mult_p(mod(i, 2) + 1), z_zero, &
416 mult_p(mod(i + 1, 2) + 1))
417 fin_p = mult_p(mod(i + 1, 2) + 1)
418 END DO
419 ELSE
420 fin_p = npq
421 END IF
422 DO i = 1, ldim
423 exp_h(1)%local_data(:, i) = real(fin_p%local_data(:, i), kind=dp)
424 exp_h(2)%local_data(:, i) = aimag(fin_p%local_data(:, i))
425 END DO
426
427 CALL cp_cfm_release(npq)
428 CALL cp_cfm_release(dpq)
429 CALL cp_cfm_release(t1)
430 CALL cp_cfm_release(t2)
431 CALL cp_cfm_release(tres)
432 CALL timestop(handle)
433
434 END SUBROUTINE exp_pade_full_complex
435
436! **************************************************************************************************
437!> \brief exponential of a complex matrix,
438!> calculated using pade approximation together with scaling and squaring
439!> \param exp_H ...
440!> \param im_part ...
441!> \param nsquare ...
442!> \param npade ...
443!> \author Florian Schiffmann (02.09)
444! **************************************************************************************************
445
446 SUBROUTINE exp_pade_only_imaginary(exp_H, im_part, nsquare, npade)
447 TYPE(cp_fm_type), DIMENSION(2) :: exp_h
448 TYPE(cp_fm_type), INTENT(IN) :: im_part
449 INTEGER, INTENT(in) :: nsquare, npade
450
451 CHARACTER(len=*), PARAMETER :: routinen = 'exp_pade_only_imaginary'
452 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
453
454 COMPLEX(KIND=dp) :: scaled, scalen
455 INTEGER :: handle, i, j, k, ldim, ndim, p, q
456 REAL(kind=dp) :: my_fac, square_fac
457 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
458 POINTER :: local_data_im
459 TYPE(cp_cfm_type) :: dpq, fin_p
460 TYPE(cp_cfm_type), DIMENSION(2) :: cmult_p
461 TYPE(cp_cfm_type), TARGET :: npq, t1
462 TYPE(cp_fm_type) :: t2, tres
463
464 CALL timeset(routinen, handle)
465 p = npade
466 q = npade !p==q seems to be necessary for the rest of the code
467
468 CALL cp_fm_get_info(im_part, local_data=local_data_im, ncol_local=ldim, nrow_global=ndim)
469 square_fac = 1.0_dp/(2.0_dp**real(nsquare, dp))
470
471 CALL cp_cfm_create(dpq, &
472 matrix_struct=im_part%matrix_struct, &
473 name="Dpq")
474
475 CALL cp_cfm_create(npq, &
476 matrix_struct=dpq%matrix_struct, &
477 name="Npq")
478
479 CALL cp_cfm_create(t1, &
480 matrix_struct=dpq%matrix_struct, &
481 name="T1")
482
483 CALL cp_fm_create(t2, &
484 matrix_struct=t1%matrix_struct, &
485 name="T2")
486
487 CALL cp_fm_create(tres, &
488 matrix_struct=t1%matrix_struct, &
489 name="Tres")
490
491! DO i=1,ldim
492! local_data_im(:,i)=local_data_im(:,i)/square_fac
493! END DO
494
495 CALL cp_fm_to_fm(im_part, t2)
496
497 CALL cp_cfm_set_all(npq, z_zero, z_one)
498 CALL cp_cfm_set_all(dpq, z_zero, z_one)
499
500 DO i = 1, ldim
501 npq%local_data(:, i) = npq%local_data(:, i) + cmplx(rzero, 0.5_dp*square_fac*local_data_im(:, i), dp)
502 dpq%local_data(:, i) = dpq%local_data(:, i) - cmplx(rzero, 0.5_dp*square_fac*local_data_im(:, i), dp)
503 END DO
504
505 IF (npade .GT. 2) THEN
506 DO j = 1, floor(npade/2.0_dp)
507 i = 2*j
508 my_fac = (-rone)**j
509 IF (i .LE. p) scalen = cmplx(my_fac*fac(p + q - i)*fac(p)/(fac(p + q)*fac(i)*fac(p - i)), 0.0_dp, dp)
510 scaled = cmplx(my_fac*fac(p + q - i)*fac(q)/(fac(p + q)*fac(i)*fac(q - i)), 0.0_dp, dp)
511 CALL parallel_gemm("N", "N", ndim, ndim, ndim, square_fac, im_part, t2, rzero, tres)
512
513 DO k = 1, ldim
514 npq%local_data(:, k) = npq%local_data(:, k) + scalen*tres%local_data(:, k)
515 dpq%local_data(:, k) = dpq%local_data(:, k) + scaled*tres%local_data(:, k)
516 END DO
517
518 IF (2*j + 1 .LE. q) THEN
519 i = 2*j + 1
520 IF (i .LE. p) scalen = cmplx(my_fac*fac(p + q - i)*fac(p)/(fac(p + q)*fac(i)*fac(p - i)), rzero, dp)
521 scaled = cmplx(-my_fac*fac(p + q - i)*fac(q)/(fac(p + q)*fac(i)*fac(q - i)), rzero, dp)
522 CALL parallel_gemm("N", "N", ndim, ndim, ndim, square_fac, im_part, tres, rzero, t2)
523
524 DO k = 1, ldim
525 npq%local_data(:, k) = npq%local_data(:, k) + scalen*cmplx(rzero, t2%local_data(:, k), dp)
526 dpq%local_data(:, k) = dpq%local_data(:, k) + scaled*cmplx(rzero, t2%local_data(:, k), dp)
527 END DO
528 END IF
529 END DO
530 END IF
531
532 CALL cp_cfm_solve(dpq, npq)
533
534 cmult_p(2) = npq
535 cmult_p(1) = t1
536 IF (nsquare .GT. 0) THEN
537 DO i = 1, nsquare
538 CALL parallel_gemm("N", "N", ndim, ndim, ndim, z_one, cmult_p(mod(i, 2) + 1), cmult_p(mod(i, 2) + 1), z_zero, &
539 cmult_p(mod(i + 1, 2) + 1))
540 fin_p = cmult_p(mod(i + 1, 2) + 1)
541 END DO
542 ELSE
543 fin_p = npq
544 END IF
545
546 DO k = 1, ldim
547 exp_h(1)%local_data(:, k) = real(fin_p%local_data(:, k), kind=dp)
548 exp_h(2)%local_data(:, k) = aimag(fin_p%local_data(:, k))
549 END DO
550
551 CALL cp_cfm_release(npq)
552 CALL cp_cfm_release(dpq)
553 CALL cp_cfm_release(t1)
554 CALL cp_fm_release(t2)
555 CALL cp_fm_release(tres)
556 CALL timestop(handle)
557 END SUBROUTINE exp_pade_only_imaginary
558
559! **************************************************************************************************
560!> \brief exponential of a real matrix,
561!> calculated using pade approximation together with scaling and squaring
562!> \param exp_H ...
563!> \param matrix ...
564!> \param nsquare ...
565!> \param npade ...
566!> \author Florian Schiffmann (02.09)
567! **************************************************************************************************
568
569 SUBROUTINE exp_pade_real(exp_H, matrix, nsquare, npade)
570 TYPE(cp_fm_type), INTENT(IN) :: exp_h, matrix
571 INTEGER, INTENT(in) :: nsquare, npade
572
573 CHARACTER(len=*), PARAMETER :: routinen = 'exp_pade_real'
574 REAL(kind=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
575
576 INTEGER :: handle, i, j, k, ldim, ndim, p, q
577 REAL(kind=dp) :: my_fac, scaled, scalen, square_fac
578 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
579 POINTER :: local_data
580 TYPE(cp_fm_type) :: dpq, fin_p
581 TYPE(cp_fm_type), DIMENSION(2) :: mult_p
582 TYPE(cp_fm_type), TARGET :: npq, t1, t2, tres
583
584 CALL timeset(routinen, handle)
585 p = npade
586 q = npade !p==q seems to be necessary for the rest of the code
587
588 CALL cp_fm_get_info(matrix, local_data=local_data, ncol_local=ldim, nrow_global=ndim)
589 square_fac = 2.0_dp**real(nsquare, dp)
590
591 CALL cp_fm_create(dpq, &
592 matrix_struct=matrix%matrix_struct, &
593 name="Dpq")
594
595 CALL cp_fm_create(npq, &
596 matrix_struct=dpq%matrix_struct, &
597 name="Npq")
598
599 CALL cp_fm_create(t1, &
600 matrix_struct=dpq%matrix_struct, &
601 name="T1")
602
603 CALL cp_fm_create(t2, &
604 matrix_struct=t1%matrix_struct, &
605 name="T2")
606
607 CALL cp_fm_create(tres, &
608 matrix_struct=t1%matrix_struct, &
609 name="Tres")
610
611 DO i = 1, ldim
612 t2%local_data(:, i) = local_data(:, i)/square_fac
613 END DO
614
615 CALL cp_fm_to_fm(t2, t1)
616 CALL cp_fm_set_all(npq, zero, one)
617 CALL cp_fm_set_all(dpq, zero, one)
618
619 DO i = 1, ldim
620 npq%local_data(:, i) = npq%local_data(:, i) + 0.5_dp*local_data(:, i)
621 dpq%local_data(:, i) = dpq%local_data(:, i) - 0.5_dp*local_data(:, i)
622 END DO
623
624 mult_p(1) = t2
625 mult_p(2) = tres
626
627 IF (npade .GE. 2) THEN
628 DO j = 2, npade
629 my_fac = (-1.0_dp)**j
630 scalen = fac(p + q - j)*fac(p)/(fac(p + q)*fac(j)*fac(p - j))
631 scaled = my_fac*fac(p + q - j)*fac(q)/(fac(p + q)*fac(j)*fac(q - j))
632 CALL parallel_gemm("N", "N", ndim, ndim, ndim, one, mult_p(mod(j, 2) + 1), t1, &
633 zero, mult_p(mod(j + 1, 2) + 1))
634
635 DO k = 1, ldim
636 npq%local_data(:, k) = npq%local_data(:, k) + scalen*mult_p(mod(j + 1, 2) + 1)%local_data(:, k)
637 dpq%local_data(:, k) = dpq%local_data(:, k) + scaled*mult_p(mod(j + 1, 2) + 1)%local_data(:, k)
638 END DO
639 END DO
640 END IF
641
642 CALL cp_fm_solve(dpq, npq)
643
644 mult_p(2) = npq
645 mult_p(1) = t1
646 IF (nsquare .GT. 0) THEN
647 DO i = 1, nsquare
648 CALL parallel_gemm("N", "N", ndim, ndim, ndim, one, mult_p(mod(i, 2) + 1), mult_p(mod(i, 2) + 1), zero, &
649 mult_p(mod(i + 1, 2) + 1))
650 END DO
651 fin_p = mult_p(mod(nsquare + 1, 2) + 1)
652 ELSE
653 fin_p = npq
654 END IF
655
656 DO k = 1, ldim
657 exp_h%local_data(:, k) = fin_p%local_data(:, k)
658 END DO
659
660 CALL cp_fm_release(npq)
661 CALL cp_fm_release(dpq)
662 CALL cp_fm_release(t1)
663 CALL cp_fm_release(t2)
664 CALL cp_fm_release(tres)
665 CALL timestop(handle)
666
667 END SUBROUTINE exp_pade_real
668
669! ***************************************************************************************************
670!> \brief exponential of a complex matrix,
671!> calculated using arnoldi subspace method (directly applies to the MOs)
672!> \param mos_old ...
673!> \param mos_new ...
674!> \param eps_exp ...
675!> \param Hre ...
676!> \param Him ...
677!> \param mos_next ...
678!> \param narn_old ...
679!> \author Florian Schiffmann (02.09)
680! **************************************************************************************************
681
682 SUBROUTINE arnoldi(mos_old, mos_new, eps_exp, Hre, Him, mos_next, narn_old)
683
684 TYPE(cp_fm_type), DIMENSION(2) :: mos_old, mos_new
685 REAL(kind=dp), INTENT(in) :: eps_exp
686 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: hre
687 TYPE(cp_fm_type), INTENT(IN) :: him
688 TYPE(cp_fm_type), DIMENSION(2), OPTIONAL :: mos_next
689 INTEGER, INTENT(inout) :: narn_old
690
691 CHARACTER(len=*), PARAMETER :: routinen = 'arnoldi'
692 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
693
694 INTEGER :: handle, i, icol_local, idim, info, j, l, &
695 mydim, nao, narnoldi, ncol_local, &
696 newdim, nmo, npade, pade_step
697 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
698 INTEGER, DIMENSION(:), POINTER :: col_indices, col_procs
699 LOGICAL :: convergence, double_col, double_row
700 REAL(dp), ALLOCATABLE, DIMENSION(:) :: last_norm, norm1, results
701 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: d, mat1, mat2, mat3, n
702 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: h_approx, h_approx_save
703 REAL(kind=dp) :: conv_norm, prefac, scaled, scalen
704 TYPE(cp_fm_struct_type), POINTER :: mo_struct, newstruct
705 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: v_mats
706 TYPE(mp_comm_type) :: col_group
707 TYPE(mp_para_env_type), POINTER :: para_env
708
709 CALL timeset(routinen, handle)
710 para_env => mos_new(1)%matrix_struct%para_env
711
712 CALL cp_fm_get_info(mos_old(1), ncol_local=ncol_local, col_indices=col_indices, &
713 nrow_global=nao, ncol_global=nmo, matrix_struct=mo_struct)
714 narnoldi = min(18, nao)
715
716 ALLOCATE (results(ncol_local))
717 ALLOCATE (norm1(ncol_local))
718 ALLOCATE (v_mats(narnoldi + 1))
719 ALLOCATE (last_norm(ncol_local))
720 ALLOCATE (h_approx(narnoldi, narnoldi, ncol_local))
721 ALLOCATE (h_approx_save(narnoldi, narnoldi, ncol_local))
722 col_procs => mo_struct%context%blacs2mpi(:, mo_struct%context%mepos(2))
723 CALL col_group%from_reordering(para_env, col_procs)
724
725 double_col = .true.
726 double_row = .false.
727 CALL cp_fm_struct_double(newstruct, mo_struct, mo_struct%context, double_col, double_row)
728 h_approx_save = rzero
729
730 DO i = 1, narnoldi + 1
731 CALL cp_fm_create(v_mats(i), matrix_struct=newstruct, &
732 name="V_mat"//cp_to_string(i))
733 END DO
734 CALL cp_fm_get_info(v_mats(1), ncol_global=newdim)
735
736 norm1 = 0.0_dp
737!$OMP PARALLEL DO PRIVATE(icol_local) DEFAULT(NONE) SHARED(V_mats,norm1,mos_old,ncol_local)
738 DO icol_local = 1, ncol_local
739 v_mats(1)%local_data(:, icol_local) = mos_old(1)%local_data(:, icol_local)
740 v_mats(1)%local_data(:, icol_local + ncol_local) = mos_old(2)%local_data(:, icol_local)
741 norm1(icol_local) = sum(v_mats(1)%local_data(:, icol_local)**2) &
742 + sum(v_mats(1)%local_data(:, icol_local + ncol_local)**2)
743 END DO
744
745 CALL col_group%sum(norm1)
746 !!! normalize the mo vectors
747 norm1(:) = sqrt(norm1(:))
748
749!$OMP PARALLEL DO PRIVATE(icol_local) DEFAULT(NONE) SHARED(V_mats,norm1,ncol_local)
750 DO icol_local = 1, ncol_local
751 v_mats(1)%local_data(:, icol_local) = v_mats(1)%local_data(:, icol_local)/norm1(icol_local)
752 v_mats(1)%local_data(:, icol_local + ncol_local) = &
753 v_mats(1)%local_data(:, icol_local + ncol_local)/norm1(icol_local)
754 END DO
755
756 ! arnoldi subspace procedure to get H_approx
757 DO i = 2, narnoldi + 1
758 !Be careful, imaginary matrix multiplied with complex. Unfortunately requires a swap of arrays afterwards
759 CALL parallel_gemm("N", "N", nao, newdim, nao, 1.0_dp, him, v_mats(i - 1), 0.0_dp, v_mats(i))
760
761!$OMP PARALLEL DO PRIVATE(icol_local) DEFAULT(NONE) SHARED(mos_new,V_mats,ncol_local,i)
762 DO icol_local = 1, ncol_local
763 mos_new(1)%local_data(:, icol_local) = v_mats(i)%local_data(:, icol_local)
764 v_mats(i)%local_data(:, icol_local) = -v_mats(i)%local_data(:, icol_local + ncol_local)
765 v_mats(i)%local_data(:, icol_local + ncol_local) = mos_new(1)%local_data(:, icol_local)
766 END DO
767
768 IF (PRESENT(hre)) THEN
769 CALL parallel_gemm("N", "N", nao, newdim, nao, 1.0_dp, hre, v_mats(i - 1), 1.0_dp, v_mats(i))
770 END IF
771
772 DO l = 1, i - 1
773!$OMP PARALLEL DO DEFAULT(NONE) SHARED(results,V_mats,ncol_local,l,i)
774 DO icol_local = 1, ncol_local
775 results(icol_local) = sum(v_mats(l)%local_data(:, icol_local)*v_mats(i)%local_data(:, icol_local)) + &
776 sum(v_mats(l)%local_data(:, icol_local + ncol_local)* &
777 v_mats(i)%local_data(:, icol_local + ncol_local))
778 END DO
779
780 CALL col_group%sum(results)
781
782!$OMP PARALLEL DO DEFAULT(NONE) SHARED(H_approx_save,V_mats,ncol_local,l,i,results)
783 DO icol_local = 1, ncol_local
784 h_approx_save(l, i - 1, icol_local) = results(icol_local)
785 v_mats(i)%local_data(:, icol_local) = v_mats(i)%local_data(:, icol_local) - &
786 results(icol_local)*v_mats(l)%local_data(:, icol_local)
787 v_mats(i)%local_data(:, icol_local + ncol_local) = &
788 v_mats(i)%local_data(:, icol_local + ncol_local) - &
789 results(icol_local)*v_mats(l)%local_data(:, icol_local + ncol_local)
790 END DO
791 END DO
792
793!$OMP PARALLEL DO DEFAULT(NONE) SHARED(ncol_local,V_mats,results,i)
794 DO icol_local = 1, ncol_local
795 results(icol_local) = sum(v_mats(i)%local_data(:, icol_local)**2) + &
796 sum(v_mats(i)%local_data(:, icol_local + ncol_local)**2)
797 END DO
798
799 CALL col_group%sum(results)
800
801 IF (i .LE. narnoldi) THEN
802
803!$OMP PARALLEL DO DEFAULT(NONE) SHARED(H_approx_save,last_norm,V_mats,ncol_local,i,results)
804 DO icol_local = 1, ncol_local
805 h_approx_save(i, i - 1, icol_local) = sqrt(results(icol_local))
806 last_norm(icol_local) = sqrt(results(icol_local))
807 v_mats(i)%local_data(:, icol_local) = v_mats(i)%local_data(:, icol_local)/sqrt(results(icol_local))
808 v_mats(i)%local_data(:, icol_local + ncol_local) = &
809 v_mats(i)%local_data(:, icol_local + ncol_local)/sqrt(results(icol_local))
810 END DO
811 ELSE
812!$OMP PARALLEL DO DEFAULT(NONE) SHARED(ncol_local,last_norm,results)
813 DO icol_local = 1, ncol_local
814 last_norm(icol_local) = sqrt(results(icol_local))
815 END DO
816 END IF
817
818 h_approx(:, :, :) = h_approx_save
819
820 ! PADE approximation for exp(H_approx), everything is done locally
821
822 convergence = .false.
823 IF (i .GE. narn_old) THEN
824 npade = 9
825 mydim = min(i, narnoldi)
826 ALLOCATE (ipivot(mydim))
827 ALLOCATE (mat1(mydim, mydim))
828 ALLOCATE (mat2(mydim, mydim))
829 ALLOCATE (mat3(mydim, mydim))
830 ALLOCATE (n(mydim, mydim))
831 ALLOCATE (d(mydim, mydim))
832 DO icol_local = 1, ncol_local
833 DO idim = 1, mydim
834 DO j = 1, mydim
835 mat1(idim, j) = h_approx(idim, j, icol_local)/16.0_dp
836 mat3(idim, j) = mat1(idim, j)
837 END DO
838 END DO
839 n = 0.0_dp
840 d = 0.0_dp
841 DO idim = 1, mydim
842 n(idim, idim) = rone
843 d(idim, idim) = rone
844 END DO
845 n(:, :) = n + 0.5_dp*mat1
846 d(:, :) = d - 0.5_dp*mat1
847 pade_step = 1
848 DO idim = 1, 4
849 pade_step = pade_step + 1
850 CALL dgemm("N", 'N', mydim, mydim, mydim, rone, mat1(1, 1), &
851 mydim, mat3(1, 1), mydim, rzero, mat2(1, 1), mydim)
852 scalen = real(fac(2*npade - pade_step)*fac(npade)/ &
853 (fac(2*npade)*fac(pade_step)*fac(npade - pade_step)), dp)
854 scaled = real((-1.0_dp)**pade_step*fac(2*npade - pade_step)*fac(npade)/ &
855 (fac(2*npade)*fac(pade_step)*fac(npade - pade_step)), dp)
856 n(:, :) = n + scalen*mat2
857 d(:, :) = d + scaled*mat2
858 pade_step = pade_step + 1
859 CALL dgemm("N", 'N', mydim, mydim, mydim, rone, mat2(1, 1), &
860 mydim, mat1(1, 1), mydim, rzero, mat3(1, 1), mydim)
861 scalen = real(fac(2*npade - pade_step)*fac(npade)/ &
862 (fac(2*npade)*fac(pade_step)*fac(npade - pade_step)), dp)
863 scaled = real((-1.0_dp)**pade_step*fac(2*npade - pade_step)*fac(npade)/ &
864 (fac(2*npade)*fac(pade_step)*fac(npade - pade_step)), dp)
865 n(:, :) = n + scalen*mat3
866 d(:, :) = d + scaled*mat3
867 END DO
868
869 CALL dgetrf(mydim, mydim, d(1, 1), mydim, ipivot, info)
870 CALL dgetrs("N", mydim, mydim, d(1, 1), mydim, ipivot, n, mydim, info)
871 CALL dgemm("N", 'N', mydim, mydim, mydim, rone, n(1, 1), mydim, n(1, 1), mydim, rzero, mat1(1, 1), mydim)
872 CALL dgemm("N", 'N', mydim, mydim, mydim, rone, mat1(1, 1), mydim, mat1(1, 1), mydim, rzero, n(1, 1), mydim)
873 CALL dgemm("N", 'N', mydim, mydim, mydim, rone, n(1, 1), mydim, n(1, 1), mydim, rzero, mat1(1, 1), mydim)
874 CALL dgemm("N", 'N', mydim, mydim, mydim, rone, mat1(1, 1), mydim, mat1(1, 1), mydim, rzero, n(1, 1), mydim)
875 DO idim = 1, mydim
876 DO j = 1, mydim
877 h_approx(idim, j, icol_local) = n(idim, j)
878 END DO
879 END DO
880 END DO
881 ! H_approx is exp(H_approx) right now, calculate new MOs and check for convergence
882 conv_norm = 0.0_dp
883 results = 0.0_dp
884 DO icol_local = 1, ncol_local
885 results(icol_local) = last_norm(icol_local)*h_approx(i - 1, 1, icol_local)
886 conv_norm = max(conv_norm, abs(results(icol_local)))
887 END DO
888
889 CALL para_env%max(conv_norm)
890
891 IF (conv_norm .LT. eps_exp .OR. i .GT. narnoldi) THEN
892
893 mos_new(1)%local_data = rzero
894 mos_new(2)%local_data = rzero
895 DO icol_local = 1, ncol_local
896 DO idim = 1, mydim
897 prefac = h_approx(idim, 1, icol_local)*norm1(icol_local)
898 mos_new(1)%local_data(:, icol_local) = mos_new(1)%local_data(:, icol_local) + &
899 v_mats(idim)%local_data(:, icol_local)*prefac
900 mos_new(2)%local_data(:, icol_local) = mos_new(2)%local_data(:, icol_local) + &
901 v_mats(idim)%local_data(:, icol_local + ncol_local)*prefac
902 END DO
903 END DO
904
905 IF (PRESENT(mos_next)) THEN
906 DO icol_local = 1, ncol_local
907 DO idim = 1, mydim
908 DO j = 1, mydim
909 n(idim, j) = h_approx(idim, j, icol_local)
910 END DO
911 END DO
912 CALL dgemm("N", 'N', mydim, mydim, mydim, rone, n(1, 1), mydim, n(1, 1), mydim, rzero, mat1(1, 1), mydim)
913 DO idim = 1, mydim
914 DO j = 1, mydim
915 h_approx(idim, j, icol_local) = mat1(idim, j)
916 END DO
917 END DO
918 END DO
919 mos_next(1)%local_data = rzero
920 mos_next(2)%local_data = rzero
921 DO icol_local = 1, ncol_local
922 DO idim = 1, mydim
923 prefac = h_approx(idim, 1, icol_local)*norm1(icol_local)
924 mos_next(1)%local_data(:, icol_local) = &
925 mos_next(1)%local_data(:, icol_local) + &
926 v_mats(idim)%local_data(:, icol_local)*prefac
927 mos_next(2)%local_data(:, icol_local) = &
928 mos_next(2)%local_data(:, icol_local) + &
929 v_mats(idim)%local_data(:, icol_local + ncol_local)*prefac
930 END DO
931 END DO
932 END IF
933 IF (conv_norm .LT. eps_exp) THEN
934 convergence = .true.
935 narn_old = i - 1
936 END IF
937 END IF
938
939 DEALLOCATE (ipivot)
940 DEALLOCATE (mat1)
941 DEALLOCATE (mat2)
942 DEALLOCATE (mat3)
943 DEALLOCATE (n)
944 DEALLOCATE (d)
945 END IF
946 IF (convergence) EXIT
947
948 END DO
949 cpwarn_if(.NOT. convergence, "ARNOLDI method did not converge")
950 !deallocate all work matrices
951
952 CALL cp_fm_release(v_mats)
953 CALL cp_fm_struct_release(newstruct)
954 CALL col_group%free()
955
956 DEALLOCATE (h_approx)
957 DEALLOCATE (h_approx_save)
958 DEALLOCATE (results)
959 DEALLOCATE (norm1)
960 DEALLOCATE (last_norm)
961 CALL timestop(handle)
962 END SUBROUTINE arnoldi
963
964END MODULE matrix_exp
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_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:55
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:82
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_scale_and_add
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
Definition cp_cfm_basic_linalg.F:240
cp_cfm_basic_linalg::cp_cfm_solve
subroutine, public cp_cfm_solve(matrix_a, general_a, determinant)
Solve the system of linear equations A*b=A_general using LU decomposition. Pay attention that both ma...
Definition cp_cfm_basic_linalg.F:732
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:121
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:155
cp_cfm_types::cp_cfm_set_all
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
Definition cp_cfm_types.F:175
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_solve
subroutine, public cp_fm_solve(matrix_a, general_a)
computes the the solution to A*b=A_general using lu decomposition
Definition cp_fm_basic_linalg.F:2338
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:167
cp_fm_basic_linalg::cp_complex_fm_gemm
subroutine, public cp_complex_fm_gemm(transa, transb, m, n, k, alpha, a_re, a_im, b_re, b_im, beta, c_re, c_im, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, c_first_row)
Convenience function. Computes the matrix multiplications needed for the multiplication of complex ma...
Definition cp_fm_basic_linalg.F:2410
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_double
subroutine, public cp_fm_struct_double(fmstruct, struct, context, col, row)
creates a struct with twice the number of blocks on each core. If matrix A has to be multiplied with ...
Definition cp_fm_struct.F:523
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1009
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:528
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:164
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
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:143
mathconstants::one
real(kind=dp), parameter, public one
Definition mathconstants.F:123
mathconstants::zero
real(kind=dp), parameter, public zero
Definition mathconstants.F:123
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
matrix_exp
Routines for calculating a complex matrix exponential.
Definition matrix_exp.F:13
matrix_exp::exp_pade_real
subroutine, public exp_pade_real(exp_h, matrix, nsquare, npade)
exponential of a real matrix, calculated using pade approximation together with scaling and squaring
Definition matrix_exp.F:570
matrix_exp::get_nsquare_norder
subroutine, public get_nsquare_norder(norm, nsquare, norder, eps_exp, method, do_emd)
optimization function for pade/taylor order and number of squaring steps
Definition matrix_exp.F:244
matrix_exp::taylor_full_complex
subroutine, public taylor_full_complex(exp_h, re_part, im_part, nsquare, ntaylor)
subroutine for general complex matrix exponentials on input a separate cp_fm_type for real and comple...
Definition matrix_exp.F:165
matrix_exp::arnoldi
subroutine, public arnoldi(mos_old, mos_new, eps_exp, hre, him, mos_next, narn_old)
exponential of a complex matrix, calculated using arnoldi subspace method (directly applies to the MO...
Definition matrix_exp.F:683
matrix_exp::exp_pade_only_imaginary
subroutine, public exp_pade_only_imaginary(exp_h, im_part, nsquare, npade)
exponential of a complex matrix, calculated using pade approximation together with scaling and squari...
Definition matrix_exp.F:447
matrix_exp::taylor_only_imaginary
subroutine, public taylor_only_imaginary(exp_h, im_matrix, nsquare, ntaylor)
specialized subroutine for purely imaginary matrix exponentials
Definition matrix_exp.F:69
matrix_exp::exp_pade_full_complex
subroutine, public exp_pade_full_complex(exp_h, re_part, im_part, nsquare, npade)
exponential of a complex matrix, calculated using pade approximation together with scaling and squari...
Definition matrix_exp.F:342
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:68
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763