(git:ccc2433)
cp_cfm_diag.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 used for collecting diagonalization schemes available for cp_cfm_type
10 !> \note
11 !> first version : only one routine right now
12 !> \author Joost VandeVondele (2003-09)
13 ! **************************************************************************************************
14 MODULE cp_cfm_diag
15  USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose, &
16  cp_cfm_gemm, &
17  cp_cfm_column_scale, &
18  cp_cfm_scale, &
19  cp_cfm_triangular_invert, &
20  cp_cfm_triangular_multiply
21  USE cp_cfm_types, ONLY: cp_cfm_get_info, &
22  cp_cfm_set_element, &
23  cp_cfm_to_cfm, &
24  cp_cfm_type
25  USE kinds, ONLY: dp
26 #if defined (__HAS_IEEE_EXCEPTIONS)
27  USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
28  ieee_set_halting_mode, &
29  ieee_all
30 #endif
31 
32 #include "../base/base_uses.f90"
33 
34  IMPLICIT NONE
35  PRIVATE
36  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_diag'
37 
38  PUBLIC :: cp_cfm_heevd, cp_cfm_geeig, cp_cfm_geeig_canon
39 
40 CONTAINS
41 
42 ! **************************************************************************************************
43 !> \brief Perform a diagonalisation of a complex matrix
44 !> \param matrix ...
45 !> \param eigenvectors ...
46 !> \param eigenvalues ...
47 !> \par History
48 !> - (De)Allocation checks updated (15.02.2011,MK)
49 !> \author Joost VandeVondele
50 ! **************************************************************************************************
51  SUBROUTINE cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
52 
53  TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
54  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
55 
56  CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_heevd'
57 
58  COMPLEX(KIND=dp), DIMENSION(:), POINTER :: work
59  COMPLEX(KIND=dp), DIMENSION(:, :), &
60  POINTER :: m
61  INTEGER :: handle, info, liwork, &
62  lrwork, lwork, n
63  INTEGER, DIMENSION(:), POINTER :: iwork
64  REAL(kind=dp), DIMENSION(:), POINTER :: rwork
65 #if defined(__SCALAPACK)
66  INTEGER, DIMENSION(9) :: descm, descv
67  COMPLEX(KIND=dp), DIMENSION(:, :), &
68  POINTER :: v
69 #if defined (__HAS_IEEE_EXCEPTIONS)
70  LOGICAL, DIMENSION(5) :: halt
71 #endif
72 #endif
73 
74  CALL timeset(routinen, handle)
75 
76  n = matrix%matrix_struct%nrow_global
77  m => matrix%local_data
78  ALLOCATE (iwork(1), rwork(1), work(1))
79  ! work space query
80  lwork = -1
81  lrwork = -1
82  liwork = -1
83 
84 #if defined(__SCALAPACK)
85  v => eigenvectors%local_data
86  descm(:) = matrix%matrix_struct%descriptor(:)
87  descv(:) = eigenvectors%matrix_struct%descriptor(:)
88  CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
89  work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
90  ! The work space query for lwork does not return always sufficiently large values.
91  ! Let's add some margin to avoid crashes.
92  lwork = ceiling(real(work(1), kind=dp)) + 1000
93  ! needed to correct for a bug in scalapack, unclear how much the right number is
94  lrwork = ceiling(real(rwork(1), kind=dp)) + 1000000
95  liwork = iwork(1)
96 #else
97  CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
98  work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
99  lwork = ceiling(real(work(1), kind=dp))
100  lrwork = ceiling(real(rwork(1), kind=dp))
101  liwork = iwork(1)
102 #endif
103 
104  DEALLOCATE (iwork, rwork, work)
105  ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
106 
107 #if defined(__SCALAPACK)
108 ! Scalapack takes advantage of IEEE754 exceptions for speedup.
109 ! Therefore, we disable floating point traps temporarily.
110 #if defined (__HAS_IEEE_EXCEPTIONS)
111  CALL ieee_get_halting_mode(ieee_all, halt)
112  CALL ieee_set_halting_mode(ieee_all, .false.)
113 #endif
114 
115  CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
116  work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
117 
118 #if defined (__HAS_IEEE_EXCEPTIONS)
119  CALL ieee_set_halting_mode(ieee_all, halt)
120 #endif
121 #else
122  CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
123  work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
124  eigenvectors%local_data = matrix%local_data
125 #endif
126 
127  DEALLOCATE (iwork, rwork, work)
128  IF (info /= 0) &
129  cpabort("Diagonalisation of a complex matrix failed")
130 
131  CALL timestop(handle)
132 
133  END SUBROUTINE cp_cfm_heevd
134 
135 ! **************************************************************************************************
136 !> \brief General Eigenvalue Problem AX = BXE
137 !> Single option version: Cholesky decomposition of B
138 !> \param amatrix ...
139 !> \param bmatrix ...
140 !> \param eigenvectors ...
141 !> \param eigenvalues ...
142 !> \param work ...
143 ! **************************************************************************************************
144  SUBROUTINE cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
145 
146  TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
147  REAL(kind=dp), DIMENSION(:) :: eigenvalues
148  TYPE(cp_cfm_type), INTENT(IN) :: work
149 
150  CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_geeig'
151 
152  INTEGER :: handle, nao, nmo
153  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
154 
155  CALL timeset(routinen, handle)
156 
157  CALL cp_cfm_get_info(amatrix, nrow_global=nao)
158  ALLOCATE (evals(nao))
159  ! Cholesky decompose S=U(T)U
160  CALL cp_cfm_cholesky_decompose(bmatrix)
161  ! Invert to get U^(-1)
162  CALL cp_cfm_triangular_invert(bmatrix)
163  ! Reduce to get U^(-T) * H * U^(-1)
164  CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
165  CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
166  ! Diagonalize
167  CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
168  ! Restore vectors C = U^(-1) * C*
169  CALL cp_cfm_triangular_multiply(bmatrix, work)
170  nmo = SIZE(eigenvalues)
171  CALL cp_cfm_to_cfm(work, eigenvectors, nmo)
172  eigenvalues(1:nmo) = evals(1:nmo)
173 
174  DEALLOCATE (evals)
175 
176  CALL timestop(handle)
177 
178  END SUBROUTINE cp_cfm_geeig
179 
180 ! **************************************************************************************************
181 !> \brief General Eigenvalue Problem AX = BXE
182 !> Use canonical orthogonalization
183 !> \param amatrix ...
184 !> \param bmatrix ...
185 !> \param eigenvectors ...
186 !> \param eigenvalues ...
187 !> \param work ...
188 !> \param epseig ...
189 ! **************************************************************************************************
190  SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
191 
192  TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
193  REAL(kind=dp), DIMENSION(:) :: eigenvalues
194  TYPE(cp_cfm_type), INTENT(IN) :: work
195  REAL(kind=dp), INTENT(IN) :: epseig
196 
197  CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_geeig_canon'
198  COMPLEX(KIND=dp), PARAMETER :: cone = cmplx(1.0_dp, 0.0_dp, kind=dp), &
199  czero = cmplx(0.0_dp, 0.0_dp, kind=dp)
200 
201  COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals
202  INTEGER :: handle, i, icol, irow, nao, nc, ncol, &
203  nmo, nx
204  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
205 
206  CALL timeset(routinen, handle)
207 
208  ! Test sizes
209  CALL cp_cfm_get_info(amatrix, nrow_global=nao)
210  nmo = SIZE(eigenvalues)
211  ALLOCATE (evals(nao), cevals(nao))
212 
213  ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
214  CALL cp_cfm_scale(-cone, bmatrix)
215  CALL cp_cfm_heevd(bmatrix, work, evals)
216  evals(:) = -evals(:)
217  nc = nao
218  DO i = 1, nao
219  IF (evals(i) < epseig) THEN
220  nc = i - 1
221  EXIT
222  END IF
223  END DO
224  cpassert(nc /= 0)
225 
226  IF (nc /= nao) THEN
227  IF (nc < nmo) THEN
228  ! Copy NULL space definition to last vectors of eigenvectors (if needed)
229  ncol = nmo - nc
230  CALL cp_cfm_to_cfm(work, eigenvectors, ncol, nc + 1, nc + 1)
231  END IF
232  ! Set NULL space in eigenvector matrix of S to zero
233  DO icol = nc + 1, nao
234  DO irow = 1, nao
235  CALL cp_cfm_set_element(work, irow, icol, czero)
236  END DO
237  END DO
238  ! Set small eigenvalues to a dummy save value
239  evals(nc + 1:nao) = 1.0_dp
240  END IF
241  ! calculate U*s**(-1/2)
242  cevals(:) = cmplx(1.0_dp/sqrt(evals(:)), 0.0_dp, kind=dp)
243  CALL cp_cfm_column_scale(work, cevals)
244  ! Reduce to get U^(-C) * H * U^(-1)
245  CALL cp_cfm_gemm("C", "N", nao, nao, nao, cone, work, amatrix, czero, bmatrix)
246  CALL cp_cfm_gemm("N", "N", nao, nao, nao, cone, bmatrix, work, czero, amatrix)
247  IF (nc /= nao) THEN
248  ! set diagonal values to save large value
249  DO icol = nc + 1, nao
250  CALL cp_cfm_set_element(amatrix, icol, icol, cmplx(10000.0_dp, 0.0_dp, kind=dp))
251  END DO
252  END IF
253  ! Diagonalize
254  CALL cp_cfm_heevd(amatrix, bmatrix, evals)
255  eigenvalues(1:nmo) = evals(1:nmo)
256  nx = min(nc, nmo)
257  ! Restore vectors C = U^(-1) * C*
258  CALL cp_cfm_gemm("N", "N", nao, nx, nc, cone, work, bmatrix, czero, eigenvectors)
259 
260  DEALLOCATE (evals)
261 
262  CALL timestop(handle)
263 
264  END SUBROUTINE cp_cfm_geeig_canon
265 
266 END MODULE cp_cfm_diag
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_gemm
subroutine, public cp_cfm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, c_first_row)
Performs one of the matrix-matrix operations: matrix_c = alpha * op1( matrix_a ) * op2( matrix_b ) + ...
Definition: cp_cfm_basic_linalg.F:567
cp_cfm_basic_linalg::cp_cfm_triangular_multiply
subroutine, public cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
Multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
Definition: cp_cfm_basic_linalg.F:1107
cp_cfm_basic_linalg::cp_cfm_cholesky_decompose
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
Definition: cp_cfm_basic_linalg.F:926
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling)
Scales columns of the full matrix by corresponding factors.
Definition: cp_cfm_basic_linalg.F:649
cp_cfm_basic_linalg::cp_cfm_triangular_invert
subroutine, public cp_cfm_triangular_invert(matrix_a, uplo, info_out)
Inverts a triangular matrix.
Definition: cp_cfm_basic_linalg.F:1189
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition: cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_geeig
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition: cp_cfm_diag.F:145
cp_cfm_diag::cp_cfm_heevd
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition: cp_cfm_diag.F:52
cp_cfm_diag::cp_cfm_geeig_canon
subroutine, public cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
General Eigenvalue Problem AX = BXE Use canonical orthogonalization.
Definition: cp_cfm_diag.F:191
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition: cp_cfm_types.F:12
cp_cfm_types::cp_cfm_set_element
subroutine, public cp_cfm_set_element(matrix, irow_global, icol_global, alpha)
Set the matrix element (irow_global,icol_global) of the full matrix to alpha.
Definition: cp_cfm_types.F:279
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition: cp_cfm_types.F:607
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34