(git:ccc2433)
matrix_exp.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Routines for calculating a complex matrix exponential.
10 !> \author Florian Schiffmann (02.09)
11 ! **************************************************************************************************
12 
13 MODULE 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 
57 CONTAINS
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  IF (.NOT. convergence) &
950  cpwarn("ARNOLDI method did not converge")
951  !deallocate all work matrices
952 
953  CALL cp_fm_release(v_mats)
954  CALL cp_fm_struct_release(newstruct)
955  CALL col_group%free()
956 
957  DEALLOCATE (h_approx)
958  DEALLOCATE (h_approx_save)
959  DEALLOCATE (results)
960  DEALLOCATE (norm1)
961  DEALLOCATE (last_norm)
962  CALL timestop(handle)
963  END SUBROUTINE arnoldi
964 
965 END 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:195
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:244
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:748
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:159
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:179
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 pay attention, both matrices are ov...
Definition: cp_fm_basic_linalg.F:2324
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:2395
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_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:536
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
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:1016
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:535
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:167
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::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::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::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::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::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::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
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
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