Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Basic linear algebra operations for complex full matrices.
10 : !> \note
11 : !> - not all functionality implemented
12 : !> \par History
13 : !> Nearly literal copy of Fawzi's routines
14 : !> \author Joost VandeVondele
15 : ! **************************************************************************************************
16 : MODULE cp_cfm_basic_linalg
17 : USE cp_blacs_env, ONLY: cp_blacs_env_type
18 : USE cp_cfm_types, ONLY: cp_cfm_create,&
19 : cp_cfm_get_info,&
20 : cp_cfm_release,&
21 : cp_cfm_to_cfm,&
22 : cp_cfm_type
23 : USE cp_fm_struct, ONLY: cp_fm_struct_equivalent
24 : USE cp_fm_types, ONLY: cp_fm_type
25 : USE cp_log_handling, ONLY: cp_to_string
26 : USE kahan_sum, ONLY: accurate_dot_product
27 : USE kinds, ONLY: dp
28 : USE mathconstants, ONLY: z_one,&
29 : z_zero
30 : USE message_passing, ONLY: mp_comm_type
31 : #include "../base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 : PRIVATE
35 :
36 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_basic_linalg'
38 :
39 : PUBLIC :: cp_cfm_column_scale, &
40 : cp_cfm_gemm, &
41 : cp_cfm_get_diag, &
42 : cp_cfm_lu_decompose, &
43 : cp_cfm_lu_invert, &
44 : cp_cfm_norm, &
45 : cp_cfm_scale, &
46 : cp_cfm_scale_and_add, &
47 : cp_cfm_scale_and_add_fm, &
48 : cp_cfm_schur_product, &
49 : cp_cfm_solve, &
50 : cp_cfm_trace, &
51 : cp_cfm_transpose, &
52 : cp_cfm_triangular_invert, &
53 : cp_cfm_triangular_multiply, &
54 : cp_cfm_rot_rows, &
55 : cp_cfm_rot_cols, &
56 : cp_cfm_det, & ! determinant of a complex matrix with correct sign
57 : cp_cfm_uplo_to_full, &
58 : cp_cfm_vectorsnorm
59 :
60 : REAL(kind=dp), EXTERNAL :: zlange, pzlange
61 :
62 : INTERFACE cp_cfm_scale
63 : MODULE PROCEDURE cp_cfm_dscale, cp_cfm_zscale
64 : END INTERFACE cp_cfm_scale
65 :
66 : ! **************************************************************************************************
67 :
68 : CONTAINS
69 :
70 : ! **************************************************************************************************
71 : !> \brief Computes the determinant (with a correct sign even in parallel environment!) of a complex square matrix
72 : !> \param matrix_a ...
73 : !> \param det_a ...
74 : !> \author A. Sinyavskiy (andrey.sinyavskiy@chem.uzh.ch)
75 : ! **************************************************************************************************
76 1500 : SUBROUTINE cp_cfm_det(matrix_a, det_a)
77 :
78 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
79 : COMPLEX(KIND=dp), INTENT(OUT) :: det_a
80 : COMPLEX(KIND=dp) :: determinant
81 : TYPE(cp_cfm_type) :: matrix_lu
82 1500 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a
83 : INTEGER :: n, i, info, P
84 1500 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
85 1500 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: diag
86 :
87 : #if defined(__parallel)
88 : INTEGER :: myprow, nprow, npcol, nrow_local, irow_local, &
89 : mypcol, ncol_local, icol_local, j
90 : INTEGER, DIMENSION(9) :: desca
91 : #endif
92 :
93 : CALL cp_cfm_create(matrix=matrix_lu, &
94 : matrix_struct=matrix_a%matrix_struct, &
95 1500 : name="A_lu"//TRIM(ADJUSTL(cp_to_string(1)))//"MATRIX")
96 1500 : CALL cp_cfm_to_cfm(matrix_a, matrix_lu)
97 :
98 1500 : a => matrix_lu%local_data
99 1500 : n = matrix_lu%matrix_struct%nrow_global
100 4500 : ALLOCATE (ipivot(n))
101 1500 : ipivot(:) = 0
102 1500 : P = 0
103 4500 : ALLOCATE (diag(n))
104 8322 : diag(:) = 0.0_dp
105 : #if defined(__parallel)
106 : ! Use LU decomposition
107 15000 : desca(:) = matrix_lu%matrix_struct%descriptor(:)
108 1500 : CALL pzgetrf(n, n, a(1, 1), 1, 1, desca, ipivot, info)
109 1500 : myprow = matrix_lu%matrix_struct%context%mepos(1)
110 1500 : mypcol = matrix_lu%matrix_struct%context%mepos(2)
111 1500 : nprow = matrix_lu%matrix_struct%context%num_pe(1)
112 1500 : npcol = matrix_lu%matrix_struct%context%num_pe(2)
113 1500 : nrow_local = matrix_lu%matrix_struct%nrow_locals(myprow)
114 1500 : ncol_local = matrix_lu%matrix_struct%ncol_locals(mypcol)
115 :
116 5031 : DO irow_local = 1, nrow_local
117 3531 : i = matrix_lu%matrix_struct%row_indices(irow_local)
118 40638 : DO icol_local = 1, ncol_local
119 35607 : j = matrix_lu%matrix_struct%col_indices(icol_local)
120 39138 : IF (i == j) diag(i) = matrix_lu%local_data(irow_local, icol_local)
121 : END DO
122 : END DO
123 15144 : CALL matrix_lu%matrix_struct%para_env%sum(diag)
124 8322 : determinant = PRODUCT(diag)
125 5031 : DO irow_local = 1, nrow_local
126 3531 : i = matrix_lu%matrix_struct%row_indices(irow_local)
127 5031 : IF (ipivot(irow_local) /= i) P = P + 1
128 : END DO
129 1500 : CALL matrix_lu%matrix_struct%para_env%sum(P)
130 : ! very important fix
131 1500 : P = P/npcol
132 : #else
133 : CALL zgetrf(n, n, a(1, 1), n, ipivot, info)
134 : DO i = 1, n
135 : diag(i) = matrix_lu%local_data(i, i)
136 : END DO
137 : determinant = PRODUCT(diag)
138 : DO i = 1, n
139 : IF (ipivot(i) /= i) P = P + 1
140 : END DO
141 : #endif
142 1500 : DEALLOCATE (ipivot)
143 1500 : DEALLOCATE (diag)
144 1500 : CALL cp_cfm_release(matrix_lu)
145 1500 : det_a = determinant*(-2*MOD(P, 2) + 1.0_dp)
146 1500 : END SUBROUTINE cp_cfm_det
147 :
148 : ! **************************************************************************************************
149 : !> \brief Computes the element-wise (Schur) product of two matrices: C = A \circ B .
150 : !> \param matrix_a the first input matrix
151 : !> \param matrix_b the second input matrix
152 : !> \param matrix_c matrix to store the result
153 : ! **************************************************************************************************
154 154 : SUBROUTINE cp_cfm_schur_product(matrix_a, matrix_b, matrix_c)
155 :
156 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, matrix_b, matrix_c
157 :
158 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_schur_product'
159 :
160 154 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b, c
161 : INTEGER :: handle, icol_local, irow_local, mypcol, &
162 : myprow, ncol_local, nrow_local
163 :
164 154 : CALL timeset(routineN, handle)
165 :
166 154 : myprow = matrix_a%matrix_struct%context%mepos(1)
167 154 : mypcol = matrix_a%matrix_struct%context%mepos(2)
168 :
169 154 : a => matrix_a%local_data
170 154 : b => matrix_b%local_data
171 154 : c => matrix_c%local_data
172 :
173 154 : nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
174 154 : ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
175 :
176 462 : DO icol_local = 1, ncol_local
177 770 : DO irow_local = 1, nrow_local
178 616 : c(irow_local, icol_local) = a(irow_local, icol_local)*b(irow_local, icol_local)
179 : END DO
180 : END DO
181 :
182 154 : CALL timestop(handle)
183 :
184 154 : END SUBROUTINE cp_cfm_schur_product
185 :
186 : ! **************************************************************************************************
187 : !> \brief Computes the element-wise (Schur) product of two matrices: C = A \circ conjg(B) .
188 : !> \param matrix_a the first input matrix
189 : !> \param matrix_b the second input matrix
190 : !> \param matrix_c matrix to store the result
191 : ! **************************************************************************************************
192 0 : SUBROUTINE cp_cfm_schur_product_cc(matrix_a, matrix_b, matrix_c)
193 :
194 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, matrix_b, matrix_c
195 :
196 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_schur_product_cc'
197 :
198 0 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b, c
199 : INTEGER :: handle, icol_local, irow_local, mypcol, &
200 : myprow, ncol_local, nrow_local
201 :
202 0 : CALL timeset(routineN, handle)
203 :
204 0 : myprow = matrix_a%matrix_struct%context%mepos(1)
205 0 : mypcol = matrix_a%matrix_struct%context%mepos(2)
206 :
207 0 : a => matrix_a%local_data
208 0 : b => matrix_b%local_data
209 0 : c => matrix_c%local_data
210 :
211 0 : nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
212 0 : ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
213 :
214 0 : DO icol_local = 1, ncol_local
215 0 : DO irow_local = 1, nrow_local
216 0 : c(irow_local, icol_local) = a(irow_local, icol_local)*CONJG(b(irow_local, icol_local))
217 : END DO
218 : END DO
219 :
220 0 : CALL timestop(handle)
221 :
222 0 : END SUBROUTINE cp_cfm_schur_product_cc
223 :
224 : ! **************************************************************************************************
225 : !> \brief Scale and add two BLACS matrices (a = alpha*a + beta*b).
226 : !> \param alpha ...
227 : !> \param matrix_a ...
228 : !> \param beta ...
229 : !> \param matrix_b ...
230 : !> \date 11.06.2001
231 : !> \author Matthias Krack
232 : !> \version 1.0
233 : !> \note
234 : !> Use explicit loops to avoid temporary arrays, as a compiler reasonably assumes that arrays
235 : !> matrix_a%local_data and matrix_b%local_data may overlap (they are referenced by pointers).
236 : !> In general case (alpha*a + beta*b) explicit loops appears to be up to two times more efficient
237 : !> than equivalent LAPACK calls (zscale, zaxpy). This is because using LAPACK calls implies
238 : !> two passes through each array, so data need to be retrieved twice if arrays are large
239 : !> enough to not fit into the processor's cache.
240 : ! **************************************************************************************************
241 659668 : SUBROUTINE cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
242 : COMPLEX(kind=dp), INTENT(in) :: alpha
243 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
244 : COMPLEX(kind=dp), INTENT(in), OPTIONAL :: beta
245 : TYPE(cp_cfm_type), INTENT(IN), OPTIONAL :: matrix_b
246 :
247 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_scale_and_add'
248 :
249 : COMPLEX(kind=dp) :: my_beta
250 659668 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b
251 : INTEGER :: handle, icol_local, irow_local, mypcol, &
252 : myprow, ncol_local, nrow_local
253 :
254 659668 : CALL timeset(routineN, handle)
255 :
256 659668 : my_beta = z_zero
257 659668 : IF (PRESENT(beta)) my_beta = beta
258 659668 : NULLIFY (a, b)
259 :
260 : ! to do: use dscal,dcopy,daxp
261 659668 : myprow = matrix_a%matrix_struct%context%mepos(1)
262 659668 : mypcol = matrix_a%matrix_struct%context%mepos(2)
263 :
264 659668 : nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
265 659668 : ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
266 :
267 659668 : a => matrix_a%local_data
268 :
269 659668 : IF (my_beta == z_zero) THEN
270 :
271 66230 : IF (alpha == z_zero) THEN
272 0 : a(:, :) = z_zero
273 66230 : ELSE IF (alpha == z_one) THEN
274 66230 : CALL timestop(handle)
275 66230 : RETURN
276 : ELSE
277 0 : a(:, :) = alpha*a(:, :)
278 : END IF
279 :
280 : ELSE
281 593438 : CPASSERT(PRESENT(matrix_b))
282 593438 : IF (matrix_a%matrix_struct%context /= matrix_b%matrix_struct%context) &
283 0 : CPABORT("matrixes must be in the same blacs context")
284 :
285 593438 : IF (cp_fm_struct_equivalent(matrix_a%matrix_struct, &
286 : matrix_b%matrix_struct)) THEN
287 :
288 593438 : b => matrix_b%local_data
289 :
290 593438 : IF (alpha == z_zero) THEN
291 20 : IF (my_beta == z_one) THEN
292 : !a(:, :) = b(:, :)
293 724 : DO icol_local = 1, ncol_local
294 49172 : DO irow_local = 1, nrow_local
295 49152 : a(irow_local, icol_local) = b(irow_local, icol_local)
296 : END DO
297 : END DO
298 : ELSE
299 : !a(:, :) = my_beta*b(:, :)
300 0 : DO icol_local = 1, ncol_local
301 0 : DO irow_local = 1, nrow_local
302 0 : a(irow_local, icol_local) = my_beta*b(irow_local, icol_local)
303 : END DO
304 : END DO
305 : END IF
306 593418 : ELSE IF (alpha == z_one) THEN
307 571875 : IF (my_beta == z_one) THEN
308 : !a(:, :) = a(:, :)+b(:, :)
309 4193997 : DO icol_local = 1, ncol_local
310 59988211 : DO irow_local = 1, nrow_local
311 59571690 : a(irow_local, icol_local) = a(irow_local, icol_local) + b(irow_local, icol_local)
312 : END DO
313 : END DO
314 : ELSE
315 : !a(:, :) = a(:, :)+my_beta*b(:, :)
316 2423339 : DO icol_local = 1, ncol_local
317 71119121 : DO irow_local = 1, nrow_local
318 70963767 : a(irow_local, icol_local) = a(irow_local, icol_local) + my_beta*b(irow_local, icol_local)
319 : END DO
320 : END DO
321 : END IF
322 : ELSE
323 : !a(:, :) = alpha*a(:, :)+my_beta*b(:, :)
324 317342 : DO icol_local = 1, ncol_local
325 8074758 : DO irow_local = 1, nrow_local
326 8053215 : a(irow_local, icol_local) = alpha*a(irow_local, icol_local) + my_beta*b(irow_local, icol_local)
327 : END DO
328 : END DO
329 : END IF
330 : ELSE
331 : CALL cp_abort(__LOCATION__, &
332 : "cp_cfm_scale_and_add is not yet implemented for cases "// &
333 0 : "where input two matrix structures are not equivalent")
334 : END IF
335 : END IF
336 593438 : CALL timestop(handle)
337 659668 : END SUBROUTINE cp_cfm_scale_and_add
338 :
339 : ! **************************************************************************************************
340 : !> \brief Scale and add two BLACS matrices (a = alpha*a + beta*b).
341 : !> where b is a real matrix (adapted from cp_cfm_scale_and_add).
342 : !> \param alpha ...
343 : !> \param matrix_a ...
344 : !> \param beta ...
345 : !> \param matrix_b ...
346 : !> \date 01.08.2014
347 : !> \author JGH
348 : !> \version 1.0
349 : ! **************************************************************************************************
350 114908 : SUBROUTINE cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
351 : COMPLEX(kind=dp), INTENT(in) :: alpha
352 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
353 : COMPLEX(kind=dp), INTENT(in) :: beta
354 : TYPE(cp_fm_type), INTENT(IN) :: matrix_b
355 :
356 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_scale_and_add_fm'
357 :
358 114908 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
359 : INTEGER :: handle, icol_local, irow_local, mypcol, &
360 : myprow, ncol_local, nrow_local
361 114908 : REAL(kind=dp), DIMENSION(:, :), POINTER :: b
362 :
363 114908 : CALL timeset(routineN, handle)
364 :
365 114908 : NULLIFY (a, b)
366 :
367 114908 : myprow = matrix_a%matrix_struct%context%mepos(1)
368 114908 : mypcol = matrix_a%matrix_struct%context%mepos(2)
369 :
370 114908 : nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
371 114908 : ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
372 :
373 114908 : a => matrix_a%local_data
374 :
375 114908 : IF (beta == z_zero) THEN
376 :
377 0 : IF (alpha == z_zero) THEN
378 0 : a(:, :) = z_zero
379 0 : ELSE IF (alpha == z_one) THEN
380 0 : CALL timestop(handle)
381 0 : RETURN
382 : ELSE
383 0 : a(:, :) = alpha*a(:, :)
384 : END IF
385 :
386 : ELSE
387 114908 : IF (matrix_a%matrix_struct%context /= matrix_b%matrix_struct%context) &
388 0 : CPABORT("matrices must be in the same blacs context")
389 :
390 114908 : IF (cp_fm_struct_equivalent(matrix_a%matrix_struct, &
391 : matrix_b%matrix_struct)) THEN
392 :
393 114908 : b => matrix_b%local_data
394 :
395 114908 : IF (alpha == z_zero) THEN
396 18849 : IF (beta == z_one) THEN
397 : !a(:, :) = b(:, :)
398 667825 : DO icol_local = 1, ncol_local
399 24712211 : DO irow_local = 1, nrow_local
400 24693398 : a(irow_local, icol_local) = b(irow_local, icol_local)
401 : END DO
402 : END DO
403 : ELSE
404 : !a(:, :) = beta*b(:, :)
405 684 : DO icol_local = 1, ncol_local
406 6516 : DO irow_local = 1, nrow_local
407 6480 : a(irow_local, icol_local) = beta*b(irow_local, icol_local)
408 : END DO
409 : END DO
410 : END IF
411 96059 : ELSE IF (alpha == z_one) THEN
412 31313 : IF (beta == z_one) THEN
413 : !a(:, :) = a(:, :)+b(:, :)
414 178296 : DO icol_local = 1, ncol_local
415 2273536 : DO irow_local = 1, nrow_local
416 2264592 : a(irow_local, icol_local) = a(irow_local, icol_local) + b(irow_local, icol_local)
417 : END DO
418 : END DO
419 : ELSE
420 : !a(:, :) = a(:, :)+beta*b(:, :)
421 757789 : DO icol_local = 1, ncol_local
422 25893447 : DO irow_local = 1, nrow_local
423 25871078 : a(irow_local, icol_local) = a(irow_local, icol_local) + beta*b(irow_local, icol_local)
424 : END DO
425 : END DO
426 : END IF
427 : ELSE
428 : !a(:, :) = alpha*a(:, :)+beta*b(:, :)
429 562290 : DO icol_local = 1, ncol_local
430 5184146 : DO irow_local = 1, nrow_local
431 5119400 : a(irow_local, icol_local) = alpha*a(irow_local, icol_local) + beta*b(irow_local, icol_local)
432 : END DO
433 : END DO
434 : END IF
435 : ELSE
436 : CALL cp_abort(__LOCATION__, &
437 : "cp_cfm_scale_and_add_fm is not yet implemented for cases "// &
438 0 : "where two input matrix structures are not equivalent")
439 : END IF
440 : END IF
441 114908 : CALL timestop(handle)
442 114908 : END SUBROUTINE cp_cfm_scale_and_add_fm
443 :
444 : ! **************************************************************************************************
445 : !> \brief Computes LU decomposition of a given matrix.
446 : !> \param matrix_a full matrix
447 : !> \param determinant determinant
448 : !> \date 11.06.2001
449 : !> \author Matthias Krack
450 : !> \version 1.0
451 : !> \note
452 : !> The actual purpose right now is to efficiently compute the determinant of a given matrix.
453 : !> The original content of the matrix is destroyed.
454 : ! **************************************************************************************************
455 0 : SUBROUTINE cp_cfm_lu_decompose(matrix_a, determinant)
456 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
457 : COMPLEX(kind=dp), INTENT(out) :: determinant
458 :
459 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_lu_decompose'
460 :
461 0 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
462 : INTEGER :: counter, handle, info, irow, nrow_global
463 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
464 :
465 : #if defined(__parallel)
466 : INTEGER :: icol, ncol_local, nrow_local
467 : INTEGER, DIMENSION(9) :: desca
468 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
469 : #else
470 : INTEGER :: lda
471 : #endif
472 :
473 0 : CALL timeset(routineN, handle)
474 :
475 0 : nrow_global = matrix_a%matrix_struct%nrow_global
476 0 : a => matrix_a%local_data
477 :
478 0 : ALLOCATE (ipivot(nrow_global))
479 : #if defined(__parallel)
480 : CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, &
481 0 : row_indices=row_indices, col_indices=col_indices)
482 :
483 0 : desca(:) = matrix_a%matrix_struct%descriptor(:)
484 0 : CALL pzgetrf(nrow_global, nrow_global, a(1, 1), 1, 1, desca, ipivot, info)
485 :
486 0 : counter = 0
487 0 : DO irow = 1, nrow_local
488 0 : IF (ipivot(irow) /= row_indices(irow)) counter = counter + 1
489 : END DO
490 :
491 0 : IF (MOD(counter, 2) == 0) THEN
492 0 : determinant = z_one
493 : ELSE
494 0 : determinant = -z_one
495 : END IF
496 :
497 : ! compute product of diagonal elements
498 : irow = 1
499 : icol = 1
500 0 : DO WHILE (irow <= nrow_local .AND. icol <= ncol_local)
501 0 : IF (row_indices(irow) < col_indices(icol)) THEN
502 0 : irow = irow + 1
503 0 : ELSE IF (row_indices(irow) > col_indices(icol)) THEN
504 0 : icol = icol + 1
505 : ELSE ! diagonal element
506 0 : determinant = determinant*a(irow, icol)
507 0 : irow = irow + 1
508 0 : icol = icol + 1
509 : END IF
510 : END DO
511 0 : CALL matrix_a%matrix_struct%para_env%prod(determinant)
512 : #else
513 : lda = SIZE(a, 1)
514 : CALL zgetrf(nrow_global, nrow_global, a(1, 1), lda, ipivot, info)
515 : counter = 0
516 : determinant = z_one
517 : DO irow = 1, nrow_global
518 : IF (ipivot(irow) /= irow) counter = counter + 1
519 : determinant = determinant*a(irow, irow)
520 : END DO
521 : IF (MOD(counter, 2) == 1) determinant = -1.0_dp*determinant
522 : #endif
523 :
524 : ! info is allowed to be zero
525 : ! this does just signal a zero diagonal element
526 0 : DEALLOCATE (ipivot)
527 :
528 0 : CALL timestop(handle)
529 0 : END SUBROUTINE cp_cfm_lu_decompose
530 :
531 : ! **************************************************************************************************
532 : !> \brief Performs one of the matrix-matrix operations:
533 : !> matrix_c = alpha * op1( matrix_a ) * op2( matrix_b ) + beta*matrix_c.
534 : !> \param transa form of op1( matrix_a ):
535 : !> op1( matrix_a ) = matrix_a, when transa == 'N' ,
536 : !> op1( matrix_a ) = matrix_a^T, when transa == 'T' ,
537 : !> op1( matrix_a ) = matrix_a^H, when transa == 'C' ,
538 : !> \param transb form of op2( matrix_b )
539 : !> \param m number of rows of the matrix op1( matrix_a )
540 : !> \param n number of columns of the matrix op2( matrix_b )
541 : !> \param k number of columns of the matrix op1( matrix_a ) as well as
542 : !> number of rows of the matrix op2( matrix_b )
543 : !> \param alpha scale factor
544 : !> \param matrix_a matrix A
545 : !> \param matrix_b matrix B
546 : !> \param beta scale factor
547 : !> \param matrix_c matrix C
548 : !> \param a_first_col (optional) the first column of the matrix_a to multiply
549 : !> \param a_first_row (optional) the first row of the matrix_a to multiply
550 : !> \param b_first_col (optional) the first column of the matrix_b to multiply
551 : !> \param b_first_row (optional) the first row of the matrix_b to multiply
552 : !> \param c_first_col (optional) the first column of the matrix_c
553 : !> \param c_first_row (optional) the first row of the matrix_c
554 : !> \date 07.06.2001
555 : !> \author Matthias Krack
556 : !> \version 1.0
557 : ! **************************************************************************************************
558 884472 : SUBROUTINE cp_cfm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, &
559 : matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, &
560 : c_first_row)
561 : CHARACTER(len=1), INTENT(IN) :: transa, transb
562 : INTEGER, INTENT(IN) :: m, n, k
563 : COMPLEX(kind=dp), INTENT(IN) :: alpha
564 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, matrix_b
565 : COMPLEX(kind=dp), INTENT(IN) :: beta
566 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_c
567 : INTEGER, INTENT(IN), OPTIONAL :: a_first_col, a_first_row, b_first_col, &
568 : b_first_row, c_first_col, c_first_row
569 :
570 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_gemm'
571 :
572 884472 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b, c
573 : INTEGER :: handle, i_a, i_b, i_c, j_a, j_b, j_c
574 : #if defined(__parallel)
575 : INTEGER, DIMENSION(9) :: desca, descb, descc
576 : #else
577 : INTEGER :: lda, ldb, ldc
578 : #endif
579 :
580 884472 : CALL timeset(routineN, handle)
581 884472 : a => matrix_a%local_data
582 884472 : b => matrix_b%local_data
583 884472 : c => matrix_c%local_data
584 :
585 884472 : i_a = 1
586 884472 : IF (PRESENT(a_first_row)) i_a = a_first_row
587 :
588 884472 : j_a = 1
589 884472 : IF (PRESENT(a_first_col)) j_a = a_first_col
590 :
591 884472 : i_b = 1
592 884472 : IF (PRESENT(b_first_row)) i_b = b_first_row
593 :
594 884472 : j_b = 1
595 884472 : IF (PRESENT(b_first_col)) j_b = b_first_col
596 :
597 884472 : i_c = 1
598 884472 : IF (PRESENT(c_first_row)) i_c = c_first_row
599 :
600 884472 : j_c = 1
601 884472 : IF (PRESENT(c_first_col)) j_c = c_first_col
602 :
603 : #if defined(__parallel)
604 8844720 : desca(:) = matrix_a%matrix_struct%descriptor(:)
605 8844720 : descb(:) = matrix_b%matrix_struct%descriptor(:)
606 8844720 : descc(:) = matrix_c%matrix_struct%descriptor(:)
607 :
608 : CALL pzgemm(transa, transb, m, n, k, alpha, a(1, 1), i_a, j_a, desca, &
609 884472 : b(1, 1), i_b, j_b, descb, beta, c(1, 1), i_c, j_c, descc)
610 : #else
611 : lda = SIZE(a, 1)
612 : ldb = SIZE(b, 1)
613 : ldc = SIZE(c, 1)
614 :
615 : ! consider zgemm3m
616 : CALL zgemm(transa, transb, m, n, k, alpha, a(i_a, j_a), &
617 : lda, b(i_b, j_b), ldb, beta, c(i_c, j_c), ldc)
618 : #endif
619 884472 : CALL timestop(handle)
620 884472 : END SUBROUTINE cp_cfm_gemm
621 :
622 : ! **************************************************************************************************
623 : !> \brief Scales columns of the full matrix by corresponding factors.
624 : !> \param matrix_a matrix to scale
625 : !> \param scaling scale factors for every column. The actual number of scaled columns is
626 : !> limited by the number of scale factors given or by the actual number of columns
627 : !> whichever is smaller.
628 : !> \author Joost VandeVondele
629 : ! **************************************************************************************************
630 53635 : SUBROUTINE cp_cfm_column_scale(matrix_a, scaling)
631 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_a
632 : COMPLEX(kind=dp), DIMENSION(:), INTENT(IN) :: scaling
633 :
634 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_column_scale'
635 :
636 53635 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
637 : INTEGER :: handle, icol_local, ncol_local, &
638 : nrow_local
639 : #if defined(__parallel)
640 53635 : INTEGER, DIMENSION(:), POINTER :: col_indices
641 : #endif
642 :
643 53635 : CALL timeset(routineN, handle)
644 :
645 53635 : a => matrix_a%local_data
646 :
647 : #if defined(__parallel)
648 53635 : CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, col_indices=col_indices)
649 53635 : ncol_local = MIN(ncol_local, SIZE(scaling))
650 :
651 1127542 : DO icol_local = 1, ncol_local
652 34128922 : a(1:nrow_local, icol_local) = scaling(col_indices(icol_local))*a(1:nrow_local, icol_local)
653 : END DO
654 : #else
655 : nrow_local = SIZE(a, 1)
656 : ncol_local = MIN(SIZE(a, 2), SIZE(scaling))
657 :
658 : DO icol_local = 1, ncol_local
659 : a(1:nrow_local, icol_local) = scaling(icol_local)*a(1:nrow_local, icol_local)
660 : END DO
661 : #endif
662 :
663 53635 : CALL timestop(handle)
664 53635 : END SUBROUTINE cp_cfm_column_scale
665 :
666 : ! **************************************************************************************************
667 : !> \brief Scales a complex matrix by a real number.
668 : !> matrix_a = alpha * matrix_b
669 : !> \param alpha scale factor
670 : !> \param matrix_a complex matrix to scale
671 : ! **************************************************************************************************
672 20771 : SUBROUTINE cp_cfm_dscale(alpha, matrix_a)
673 : REAL(kind=dp), INTENT(IN) :: alpha
674 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_a
675 :
676 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_dscale'
677 :
678 20771 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
679 : INTEGER :: handle
680 :
681 20771 : CALL timeset(routineN, handle)
682 :
683 20771 : NULLIFY (a)
684 :
685 20771 : a => matrix_a%local_data
686 :
687 62313 : CALL zdscal(SIZE(a), alpha, a(1, 1), 1)
688 :
689 20771 : CALL timestop(handle)
690 20771 : END SUBROUTINE cp_cfm_dscale
691 :
692 : ! **************************************************************************************************
693 : !> \brief Scales a complex matrix by a complex number.
694 : !> matrix_a = alpha * matrix_b
695 : !> \param alpha scale factor
696 : !> \param matrix_a complex matrix to scale
697 : !> \note
698 : !> use cp_fm_set_all to zero (avoids problems with nan)
699 : ! **************************************************************************************************
700 38755 : SUBROUTINE cp_cfm_zscale(alpha, matrix_a)
701 : COMPLEX(kind=dp), INTENT(IN) :: alpha
702 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_a
703 :
704 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_zscale'
705 :
706 38755 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
707 : INTEGER :: handle
708 :
709 38755 : CALL timeset(routineN, handle)
710 :
711 38755 : NULLIFY (a)
712 :
713 38755 : a => matrix_a%local_data
714 :
715 6370108 : a(:, :) = alpha*a(:, :)
716 :
717 38755 : CALL timestop(handle)
718 38755 : END SUBROUTINE cp_cfm_zscale
719 :
720 : ! **************************************************************************************************
721 : !> \brief Solve the system of linear equations A*b=A_general using LU decomposition.
722 : !> Pay attention that both matrices are overwritten on exit and that
723 : !> the result is stored into the matrix 'general_a'.
724 : !> \param matrix_a matrix A (overwritten on exit)
725 : !> \param general_a (input) matrix A_general, (output) matrix B
726 : !> \param determinant (optional) determinant
727 : !> \author Florian Schiffmann
728 : ! **************************************************************************************************
729 7146 : SUBROUTINE cp_cfm_solve(matrix_a, general_a, determinant)
730 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, general_a
731 : COMPLEX(kind=dp), OPTIONAL :: determinant
732 :
733 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_solve'
734 :
735 7146 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, a_general
736 : INTEGER :: counter, handle, info, irow, nrow_global
737 7146 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
738 :
739 : #if defined(__parallel)
740 : INTEGER :: icol, ncol_local, nrow_local
741 : INTEGER, DIMENSION(9) :: desca, descb
742 7146 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
743 : #else
744 : INTEGER :: lda, ldb
745 : #endif
746 :
747 7146 : CALL timeset(routineN, handle)
748 :
749 7146 : a => matrix_a%local_data
750 7146 : a_general => general_a%local_data
751 7146 : nrow_global = matrix_a%matrix_struct%nrow_global
752 21438 : ALLOCATE (ipivot(nrow_global))
753 :
754 : #if defined(__parallel)
755 71460 : desca(:) = matrix_a%matrix_struct%descriptor(:)
756 71460 : descb(:) = general_a%matrix_struct%descriptor(:)
757 7146 : CALL pzgetrf(nrow_global, nrow_global, a(1, 1), 1, 1, desca, ipivot, info)
758 7146 : IF (PRESENT(determinant)) THEN
759 : CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, &
760 6418 : row_indices=row_indices, col_indices=col_indices)
761 :
762 6418 : counter = 0
763 19302 : DO irow = 1, nrow_local
764 19302 : IF (ipivot(irow) /= row_indices(irow)) counter = counter + 1
765 : END DO
766 :
767 6418 : IF (MOD(counter, 2) == 0) THEN
768 6408 : determinant = z_one
769 : ELSE
770 10 : determinant = -z_one
771 : END IF
772 :
773 : ! compute product of diagonal elements
774 : irow = 1
775 : icol = 1
776 28941 : DO WHILE (irow <= nrow_local .AND. icol <= ncol_local)
777 28941 : IF (row_indices(irow) < col_indices(icol)) THEN
778 0 : irow = irow + 1
779 22523 : ELSE IF (row_indices(irow) > col_indices(icol)) THEN
780 9639 : icol = icol + 1
781 : ELSE ! diagonal element
782 12884 : determinant = determinant*a(irow, icol)
783 12884 : irow = irow + 1
784 12884 : icol = icol + 1
785 : END IF
786 : END DO
787 6418 : CALL matrix_a%matrix_struct%para_env%prod(determinant)
788 : END IF
789 :
790 : CALL pzgetrs("N", nrow_global, nrow_global, a(1, 1), 1, 1, desca, &
791 7146 : ipivot, a_general(1, 1), 1, 1, descb, info)
792 : #else
793 : lda = SIZE(a, 1)
794 : ldb = SIZE(a_general, 1)
795 : CALL zgetrf(nrow_global, nrow_global, a(1, 1), lda, ipivot, info)
796 : IF (PRESENT(determinant)) THEN
797 : counter = 0
798 : determinant = z_one
799 : DO irow = 1, nrow_global
800 : IF (ipivot(irow) /= irow) counter = counter + 1
801 : determinant = determinant*a(irow, irow)
802 : END DO
803 : IF (MOD(counter, 2) == 1) determinant = -1.0_dp*determinant
804 : END IF
805 : CALL zgetrs("N", nrow_global, nrow_global, a(1, 1), lda, ipivot, a_general(1, 1), ldb, info)
806 : #endif
807 :
808 : ! info is allowed to be zero
809 : ! this does just signal a zero diagonal element
810 7146 : DEALLOCATE (ipivot)
811 7146 : CALL timestop(handle)
812 :
813 7146 : END SUBROUTINE cp_cfm_solve
814 :
815 : ! **************************************************************************************************
816 : !> \brief Inverts a matrix using LU decomposition. The input matrix will be overwritten.
817 : !> \param matrix input a general square non-singular matrix, outputs its inverse
818 : !> \param info_out optional, if present outputs the info from (p)zgetri
819 : !> \author Lianheng Tong
820 : ! **************************************************************************************************
821 120050 : SUBROUTINE cp_cfm_lu_invert(matrix, info_out)
822 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
823 : INTEGER, INTENT(out), OPTIONAL :: info_out
824 :
825 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_lu_invert'
826 :
827 120050 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
828 : COMPLEX(kind=dp), DIMENSION(1) :: work1
829 120050 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: mat
830 : INTEGER :: handle, info, lwork, nrows_global
831 120050 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
832 :
833 : #if defined(__parallel)
834 : INTEGER :: liwork
835 120050 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
836 : INTEGER, DIMENSION(1) :: iwork1
837 : INTEGER, DIMENSION(9) :: desca
838 : #else
839 : INTEGER :: lda
840 : #endif
841 :
842 120050 : CALL timeset(routineN, handle)
843 :
844 120050 : mat => matrix%local_data
845 120050 : nrows_global = matrix%matrix_struct%nrow_global
846 120050 : CPASSERT(nrows_global == matrix%matrix_struct%ncol_global)
847 360150 : ALLOCATE (ipivot(nrows_global))
848 :
849 : ! do LU decomposition
850 : #if defined(__parallel)
851 1200500 : desca = matrix%matrix_struct%descriptor
852 : CALL pzgetrf(nrows_global, nrows_global, &
853 120050 : mat(1, 1), 1, 1, desca, ipivot, info)
854 : #else
855 : lda = SIZE(mat, 1)
856 : CALL zgetrf(nrows_global, nrows_global, &
857 : mat(1, 1), lda, ipivot, info)
858 : #endif
859 120050 : IF (info /= 0) THEN
860 0 : CALL cp_abort(__LOCATION__, "LU decomposition has failed")
861 : END IF
862 :
863 : ! do inversion
864 : #if defined(__parallel)
865 : CALL pzgetri(nrows_global, mat(1, 1), 1, 1, desca, &
866 120050 : ipivot, work1, -1, iwork1, -1, info)
867 120050 : lwork = INT(work1(1))
868 120050 : liwork = INT(iwork1(1))
869 360150 : ALLOCATE (work(lwork))
870 360150 : ALLOCATE (iwork(liwork))
871 : CALL pzgetri(nrows_global, mat(1, 1), 1, 1, desca, &
872 120050 : ipivot, work, lwork, iwork, liwork, info)
873 120050 : DEALLOCATE (iwork)
874 : #else
875 : CALL zgetri(nrows_global, mat(1, 1), lda, ipivot, work1, -1, info)
876 : lwork = INT(work1(1))
877 : ALLOCATE (work(lwork))
878 : CALL zgetri(nrows_global, mat(1, 1), lda, ipivot, work, lwork, info)
879 : #endif
880 120050 : DEALLOCATE (work)
881 120050 : DEALLOCATE (ipivot)
882 :
883 120050 : IF (PRESENT(info_out)) THEN
884 0 : info_out = info
885 : ELSE
886 120050 : IF (info /= 0) &
887 0 : CALL cp_abort(__LOCATION__, "LU inversion has failed")
888 : END IF
889 :
890 120050 : CALL timestop(handle)
891 :
892 120050 : END SUBROUTINE cp_cfm_lu_invert
893 :
894 : ! **************************************************************************************************
895 : !> \brief Returns the trace of matrix_a^T matrix_b, i.e
896 : !> sum_{i,j}(matrix_a(i,j)*matrix_b(i,j)) .
897 : !> \param matrix_a a complex matrix
898 : !> \param matrix_b another complex matrix
899 : !> \param trace value of the trace operator
900 : !> \par History
901 : !> * 09.2017 created [Sergey Chulkov]
902 : !> \author Sergey Chulkov
903 : !> \note
904 : !> Based on the subroutine cp_fm_trace(). Note the transposition of matrix_a!
905 : ! **************************************************************************************************
906 166414 : SUBROUTINE cp_cfm_trace(matrix_a, matrix_b, trace)
907 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, matrix_b
908 : COMPLEX(kind=dp), INTENT(out) :: trace
909 :
910 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_trace'
911 :
912 : INTEGER :: handle, mypcol, myprow, ncol_local, &
913 : npcol, nprow, nrow_local
914 : TYPE(cp_blacs_env_type), POINTER :: context
915 : TYPE(mp_comm_type) :: group
916 :
917 166414 : CALL timeset(routineN, handle)
918 :
919 166414 : context => matrix_a%matrix_struct%context
920 166414 : myprow = context%mepos(1)
921 166414 : mypcol = context%mepos(2)
922 166414 : nprow = context%num_pe(1)
923 166414 : npcol = context%num_pe(2)
924 :
925 166414 : group = matrix_a%matrix_struct%para_env
926 :
927 166414 : nrow_local = MIN(matrix_a%matrix_struct%nrow_locals(myprow), matrix_b%matrix_struct%nrow_locals(myprow))
928 166414 : ncol_local = MIN(matrix_a%matrix_struct%ncol_locals(mypcol), matrix_b%matrix_struct%ncol_locals(mypcol))
929 :
930 : ! compute an accurate dot-product
931 : trace = accurate_dot_product(matrix_a%local_data(1:nrow_local, 1:ncol_local), &
932 166414 : matrix_b%local_data(1:nrow_local, 1:ncol_local))
933 :
934 166414 : CALL group%sum(trace)
935 :
936 166414 : CALL timestop(handle)
937 :
938 166414 : END SUBROUTINE cp_cfm_trace
939 :
940 : ! **************************************************************************************************
941 : !> \brief Multiplies in place by a triangular matrix:
942 : !> matrix_b = alpha op(triangular_matrix) matrix_b
943 : !> or (if side='R')
944 : !> matrix_b = alpha matrix_b op(triangular_matrix)
945 : !> op(triangular_matrix) is:
946 : !> triangular_matrix (if transa="N" and invert_tr=.false.)
947 : !> triangular_matrix^T (if transa="T" and invert_tr=.false.)
948 : !> triangular_matrix^H (if transa="C" and invert_tr=.false.)
949 : !> triangular_matrix^(-1) (if transa="N" and invert_tr=.true.)
950 : !> triangular_matrix^(-T) (if transa="T" and invert_tr=.true.)
951 : !> triangular_matrix^(-H) (if transa="C" and invert_tr=.true.)
952 : !> \param triangular_matrix the triangular matrix that multiplies the other
953 : !> \param matrix_b the matrix that gets multiplied and stores the result
954 : !> \param side on which side of matrix_b stays op(triangular_matrix)
955 : !> (defaults to 'L')
956 : !> \param transa_tr ...
957 : !> \param invert_tr if the triangular matrix should be inverted
958 : !> (defaults to false)
959 : !> \param uplo_tr if triangular_matrix is stored in the upper ('U') or
960 : !> lower ('L') triangle (defaults to 'U')
961 : !> \param unit_diag_tr if the diagonal elements of triangular_matrix should
962 : !> be assumed to be 1 (defaults to false)
963 : !> \param n_rows the number of rows of the result (defaults to
964 : !> size(matrix_b,1))
965 : !> \param n_cols the number of columns of the result (defaults to
966 : !> size(matrix_b,2))
967 : !> \param alpha ...
968 : !> \par History
969 : !> 08.2002 created [fawzi]
970 : !> \author Fawzi Mohamed
971 : !> \note
972 : !> needs an mpi env
973 : ! **************************************************************************************************
974 224326 : SUBROUTINE cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, &
975 : transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, &
976 : alpha)
977 : TYPE(cp_cfm_type), INTENT(IN) :: triangular_matrix, matrix_b
978 : CHARACTER, INTENT(in), OPTIONAL :: side, transa_tr
979 : LOGICAL, INTENT(in), OPTIONAL :: invert_tr
980 : CHARACTER, INTENT(in), OPTIONAL :: uplo_tr
981 : LOGICAL, INTENT(in), OPTIONAL :: unit_diag_tr
982 : INTEGER, INTENT(in), OPTIONAL :: n_rows, n_cols
983 : COMPLEX(kind=dp), INTENT(in), OPTIONAL :: alpha
984 :
985 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_triangular_multiply'
986 :
987 : CHARACTER :: side_char, transa, unit_diag, uplo
988 : COMPLEX(kind=dp) :: al
989 : INTEGER :: handle, m, n
990 : LOGICAL :: invert
991 :
992 112163 : CALL timeset(routineN, handle)
993 112163 : side_char = 'L'
994 112163 : unit_diag = 'N'
995 112163 : uplo = 'U'
996 112163 : transa = 'N'
997 112163 : invert = .FALSE.
998 112163 : al = z_one
999 112163 : CALL cp_cfm_get_info(matrix_b, nrow_global=m, ncol_global=n)
1000 112163 : IF (PRESENT(side)) side_char = side
1001 112163 : IF (PRESENT(invert_tr)) invert = invert_tr
1002 112163 : IF (PRESENT(uplo_tr)) uplo = uplo_tr
1003 112163 : IF (PRESENT(unit_diag_tr)) THEN
1004 0 : IF (unit_diag_tr) THEN
1005 0 : unit_diag = 'U'
1006 : ELSE
1007 : unit_diag = 'N'
1008 : END IF
1009 : END IF
1010 112163 : IF (PRESENT(transa_tr)) transa = transa_tr
1011 112163 : IF (PRESENT(alpha)) al = alpha
1012 112163 : IF (PRESENT(n_rows)) m = n_rows
1013 112163 : IF (PRESENT(n_cols)) n = n_cols
1014 :
1015 112163 : IF (invert) THEN
1016 :
1017 : #if defined(__parallel)
1018 : CALL pztrsm(side_char, uplo, transa, unit_diag, m, n, al, &
1019 : triangular_matrix%local_data(1, 1), 1, 1, &
1020 : triangular_matrix%matrix_struct%descriptor, &
1021 : matrix_b%local_data(1, 1), 1, 1, &
1022 2129 : matrix_b%matrix_struct%descriptor(1))
1023 : #else
1024 : CALL ztrsm(side_char, uplo, transa, unit_diag, m, n, al, &
1025 : triangular_matrix%local_data(1, 1), &
1026 : SIZE(triangular_matrix%local_data, 1), &
1027 : matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1))
1028 : #endif
1029 :
1030 : ELSE
1031 :
1032 : #if defined(__parallel)
1033 : CALL pztrmm(side_char, uplo, transa, unit_diag, m, n, al, &
1034 : triangular_matrix%local_data(1, 1), 1, 1, &
1035 : triangular_matrix%matrix_struct%descriptor, &
1036 : matrix_b%local_data(1, 1), 1, 1, &
1037 110034 : matrix_b%matrix_struct%descriptor(1))
1038 : #else
1039 : CALL ztrmm(side_char, uplo, transa, unit_diag, m, n, al, &
1040 : triangular_matrix%local_data(1, 1), &
1041 : SIZE(triangular_matrix%local_data, 1), &
1042 : matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1))
1043 : #endif
1044 :
1045 : END IF
1046 :
1047 112163 : CALL timestop(handle)
1048 :
1049 112163 : END SUBROUTINE cp_cfm_triangular_multiply
1050 :
1051 : ! **************************************************************************************************
1052 : !> \brief Inverts a triangular matrix.
1053 : !> \param matrix_a ...
1054 : !> \param uplo ...
1055 : !> \param info_out ...
1056 : !> \author MI
1057 : ! **************************************************************************************************
1058 36678 : SUBROUTINE cp_cfm_triangular_invert(matrix_a, uplo, info_out)
1059 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
1060 : CHARACTER, INTENT(in), OPTIONAL :: uplo
1061 : INTEGER, INTENT(out), OPTIONAL :: info_out
1062 :
1063 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_triangular_invert'
1064 :
1065 : CHARACTER :: unit_diag, my_uplo
1066 : INTEGER :: handle, info, ncol_global
1067 : COMPLEX(kind=dp), DIMENSION(:, :), &
1068 36678 : POINTER :: a
1069 : #if defined(__parallel)
1070 : INTEGER, DIMENSION(9) :: desca
1071 : #endif
1072 :
1073 36678 : CALL timeset(routineN, handle)
1074 :
1075 36678 : unit_diag = 'N'
1076 36678 : my_uplo = 'U'
1077 36678 : IF (PRESENT(uplo)) my_uplo = uplo
1078 :
1079 36678 : ncol_global = matrix_a%matrix_struct%ncol_global
1080 :
1081 36678 : a => matrix_a%local_data
1082 :
1083 : #if defined(__parallel)
1084 366780 : desca(:) = matrix_a%matrix_struct%descriptor(:)
1085 36678 : CALL pztrtri(my_uplo, unit_diag, ncol_global, a(1, 1), 1, 1, desca, info)
1086 : #else
1087 : CALL ztrtri(my_uplo, unit_diag, ncol_global, a(1, 1), ncol_global, info)
1088 : #endif
1089 :
1090 36678 : IF (PRESENT(info_out)) THEN
1091 0 : info_out = info
1092 : ELSE
1093 36678 : IF (info /= 0) &
1094 : CALL cp_abort(__LOCATION__, &
1095 0 : "triangular invert failed: matrix is not positive definite or ill-conditioned")
1096 : END IF
1097 :
1098 36678 : CALL timestop(handle)
1099 36678 : END SUBROUTINE cp_cfm_triangular_invert
1100 :
1101 : ! **************************************************************************************************
1102 : !> \brief Transposes a BLACS distributed complex matrix.
1103 : !> \param matrix input matrix
1104 : !> \param trans 'T' for transpose, 'C' for Hermitian conjugate
1105 : !> \param matrixt output matrix
1106 : !> \author Lianheng Tong
1107 : ! **************************************************************************************************
1108 30981 : SUBROUTINE cp_cfm_transpose(matrix, trans, matrixt)
1109 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1110 : CHARACTER, INTENT(in) :: trans
1111 : TYPE(cp_cfm_type), INTENT(IN) :: matrixt
1112 :
1113 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_transpose'
1114 :
1115 30981 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa, cc
1116 : INTEGER :: handle, ncol_global, nrow_global
1117 : #if defined(__parallel)
1118 : INTEGER, DIMENSION(9) :: desca, descc
1119 : #elif !defined(__MKL)
1120 : INTEGER :: ii, jj
1121 : #endif
1122 :
1123 30981 : CALL timeset(routineN, handle)
1124 :
1125 30981 : nrow_global = matrix%matrix_struct%nrow_global
1126 30981 : ncol_global = matrix%matrix_struct%ncol_global
1127 :
1128 30981 : CPASSERT(matrixt%matrix_struct%nrow_global == ncol_global)
1129 30981 : CPASSERT(matrixt%matrix_struct%ncol_global == nrow_global)
1130 :
1131 30981 : aa => matrix%local_data
1132 30981 : cc => matrixt%local_data
1133 :
1134 : #if defined(__parallel)
1135 309810 : desca = matrix%matrix_struct%descriptor
1136 309810 : descc = matrixt%matrix_struct%descriptor
1137 13338 : SELECT CASE (trans)
1138 : CASE ('T')
1139 : CALL pztranu(nrow_global, ncol_global, &
1140 : z_one, aa(1, 1), 1, 1, desca, &
1141 13338 : z_zero, cc(1, 1), 1, 1, descc)
1142 : CASE ('C')
1143 : CALL pztranc(nrow_global, ncol_global, &
1144 : z_one, aa(1, 1), 1, 1, desca, &
1145 17643 : z_zero, cc(1, 1), 1, 1, descc)
1146 : CASE DEFAULT
1147 30981 : CPABORT("trans only accepts 'T' or 'C'")
1148 : END SELECT
1149 : #elif defined(__MKL)
1150 : CALL mkl_zomatcopy('C', trans, nrow_global, ncol_global, 1.0_dp, aa(1, 1), nrow_global, cc(1, 1), ncol_global)
1151 : #else
1152 : SELECT CASE (trans)
1153 : CASE ('T')
1154 : DO jj = 1, ncol_global
1155 : DO ii = 1, nrow_global
1156 : cc(ii, jj) = aa(jj, ii)
1157 : END DO
1158 : END DO
1159 : CASE ('C')
1160 : DO jj = 1, ncol_global
1161 : DO ii = 1, nrow_global
1162 : cc(ii, jj) = CONJG(aa(jj, ii))
1163 : END DO
1164 : END DO
1165 : CASE DEFAULT
1166 : CPABORT("trans only accepts 'T' or 'C'")
1167 : END SELECT
1168 : #endif
1169 :
1170 30981 : CALL timestop(handle)
1171 30981 : END SUBROUTINE cp_cfm_transpose
1172 :
1173 : ! **************************************************************************************************
1174 : !> \brief Norm of matrix using (p)zlange.
1175 : !> \param matrix input a general matrix
1176 : !> \param mode 'M' max abs element value,
1177 : !> '1' or 'O' one norm, i.e. maximum column sum,
1178 : !> 'I' infinity norm, i.e. maximum row sum,
1179 : !> 'F' or 'E' Frobenius norm, i.e. sqrt of sum of all squares of elements
1180 : !> \return the norm according to mode
1181 : !> \author Lianheng Tong
1182 : ! **************************************************************************************************
1183 237270 : FUNCTION cp_cfm_norm(matrix, mode) RESULT(res)
1184 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1185 : CHARACTER, INTENT(IN) :: mode
1186 : REAL(kind=dp) :: res
1187 :
1188 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_norm'
1189 :
1190 237270 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa
1191 : INTEGER :: handle, lwork, ncols, ncols_local, &
1192 : nrows, nrows_local
1193 237270 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
1194 :
1195 : #if defined(__parallel)
1196 : INTEGER, DIMENSION(9) :: desca
1197 : #else
1198 : INTEGER :: lda
1199 : #endif
1200 :
1201 237270 : CALL timeset(routineN, handle)
1202 :
1203 : CALL cp_cfm_get_info(matrix=matrix, &
1204 : nrow_global=nrows, &
1205 : ncol_global=ncols, &
1206 : nrow_local=nrows_local, &
1207 237270 : ncol_local=ncols_local)
1208 237270 : aa => matrix%local_data
1209 :
1210 : SELECT CASE (mode)
1211 : CASE ('M', 'm')
1212 0 : lwork = 1
1213 : CASE ('1', 'O', 'o')
1214 : #if defined(__parallel)
1215 0 : lwork = ncols_local
1216 : #else
1217 : lwork = 1
1218 : #endif
1219 : CASE ('I', 'i')
1220 : #if defined(__parallel)
1221 12 : lwork = nrows_local
1222 : #else
1223 : lwork = nrows
1224 : #endif
1225 : CASE ('F', 'f', 'E', 'e')
1226 0 : lwork = 1
1227 : CASE DEFAULT
1228 237270 : CPABORT("mode input is not valid")
1229 : END SELECT
1230 :
1231 711810 : ALLOCATE (work(lwork))
1232 :
1233 : #if defined(__parallel)
1234 2372700 : desca = matrix%matrix_struct%descriptor
1235 237270 : res = pzlange(mode, nrows, ncols, aa(1, 1), 1, 1, desca, work)
1236 : #else
1237 : lda = SIZE(aa, 1)
1238 : res = zlange(mode, nrows, ncols, aa(1, 1), lda, work)
1239 : #endif
1240 :
1241 237270 : DEALLOCATE (work)
1242 237270 : CALL timestop(handle)
1243 237270 : END FUNCTION cp_cfm_norm
1244 :
1245 : ! **************************************************************************************************
1246 : !> \brief Applies a planar rotation defined by cs and sn to the i'th and j'th rows.
1247 : !> \param matrix ...
1248 : !> \param irow ...
1249 : !> \param jrow ...
1250 : !> \param cs cosine of the rotation angle
1251 : !> \param sn sinus of the rotation angle
1252 : !> \author Ole Schuett
1253 : ! **************************************************************************************************
1254 375120 : SUBROUTINE cp_cfm_rot_rows(matrix, irow, jrow, cs, sn)
1255 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1256 : INTEGER, INTENT(IN) :: irow, jrow
1257 : REAL(dp), INTENT(IN) :: cs, sn
1258 :
1259 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_rot_rows'
1260 : INTEGER :: handle, ncol
1261 : COMPLEX(KIND=dp) :: sn_cmplx
1262 :
1263 : #if defined(__parallel)
1264 : INTEGER :: info, lwork
1265 : INTEGER, DIMENSION(9) :: desc
1266 375120 : REAL(dp), DIMENSION(:), ALLOCATABLE :: work
1267 : #endif
1268 375120 : CALL timeset(routineN, handle)
1269 375120 : CALL cp_cfm_get_info(matrix, ncol_global=ncol)
1270 375120 : sn_cmplx = CMPLX(sn, 0.0_dp, dp)
1271 : #if defined(__parallel)
1272 375120 : IF (1 /= matrix%matrix_struct%context%n_pid) THEN
1273 375120 : lwork = 2*ncol + 1
1274 1125360 : ALLOCATE (work(lwork))
1275 3751200 : desc(:) = matrix%matrix_struct%descriptor(:)
1276 375120 : info = 0
1277 : CALL pzrot(ncol, &
1278 : matrix%local_data(1, 1), irow, 1, desc, ncol, &
1279 : matrix%local_data(1, 1), jrow, 1, desc, ncol, &
1280 375120 : cs, sn_cmplx, work, lwork, info)
1281 375120 : CPASSERT(info == 0)
1282 375120 : DEALLOCATE (work)
1283 : ELSE
1284 : #endif
1285 0 : CALL zrot(ncol, matrix%local_data(irow, 1), ncol, matrix%local_data(jrow, 1), ncol, cs, sn_cmplx)
1286 : #if defined(__parallel)
1287 : END IF
1288 : #endif
1289 375120 : CALL timestop(handle)
1290 375120 : END SUBROUTINE cp_cfm_rot_rows
1291 :
1292 : ! **************************************************************************************************
1293 : !> \brief Applies a planar rotation defined by cs and sn to the i'th and j'th columnns.
1294 : !> \param matrix ...
1295 : !> \param icol ...
1296 : !> \param jcol ...
1297 : !> \param cs cosine of the rotation angle
1298 : !> \param sn sinus of the rotation angle
1299 : !> \author Ole Schuett
1300 : ! **************************************************************************************************
1301 422760 : SUBROUTINE cp_cfm_rot_cols(matrix, icol, jcol, cs, sn)
1302 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1303 : INTEGER, INTENT(IN) :: icol, jcol
1304 : REAL(dp), INTENT(IN) :: cs, sn
1305 :
1306 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_rot_cols'
1307 : INTEGER :: handle, nrow
1308 : COMPLEX(KIND=dp) :: sn_cmplx
1309 :
1310 : #if defined(__parallel)
1311 : INTEGER :: info, lwork
1312 : INTEGER, DIMENSION(9) :: desc
1313 422760 : REAL(dp), DIMENSION(:), ALLOCATABLE :: work
1314 : #endif
1315 422760 : CALL timeset(routineN, handle)
1316 422760 : CALL cp_cfm_get_info(matrix, nrow_global=nrow)
1317 422760 : sn_cmplx = CMPLX(sn, 0.0_dp, dp)
1318 : #if defined(__parallel)
1319 422760 : IF (1 /= matrix%matrix_struct%context%n_pid) THEN
1320 422760 : lwork = 2*nrow + 1
1321 1268280 : ALLOCATE (work(lwork))
1322 4227600 : desc(:) = matrix%matrix_struct%descriptor(:)
1323 422760 : info = 0
1324 : CALL pzrot(nrow, &
1325 : matrix%local_data(1, 1), 1, icol, desc, 1, &
1326 : matrix%local_data(1, 1), 1, jcol, desc, 1, &
1327 422760 : cs, sn_cmplx, work, lwork, info)
1328 422760 : CPASSERT(info == 0)
1329 422760 : DEALLOCATE (work)
1330 : ELSE
1331 : #endif
1332 0 : CALL zrot(nrow, matrix%local_data(1, icol), 1, matrix%local_data(1, jcol), 1, cs, sn_cmplx)
1333 : #if defined(__parallel)
1334 : END IF
1335 : #endif
1336 422760 : CALL timestop(handle)
1337 422760 : END SUBROUTINE cp_cfm_rot_cols
1338 :
1339 : ! **************************************************************************************************
1340 : !> \brief ...
1341 : !> \param matrix ...
1342 : !> \param workspace ...
1343 : !> \param uplo triangular format; defaults to 'U'
1344 : !> \par History
1345 : !> 12.2024 Added optional workspace as input [Rocco Meli]
1346 : !> \author Jan Wilhelm
1347 : ! **************************************************************************************************
1348 11906 : SUBROUTINE cp_cfm_uplo_to_full(matrix, workspace, uplo)
1349 :
1350 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1351 : TYPE(cp_cfm_type), INTENT(IN), OPTIONAL :: workspace
1352 : CHARACTER, INTENT(IN), OPTIONAL :: uplo
1353 :
1354 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_cfm_uplo_to_full'
1355 :
1356 : CHARACTER :: myuplo
1357 : INTEGER :: handle, i_global, iiB, j_global, jjB, &
1358 : ncol_local, nrow_local
1359 5953 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1360 : TYPE(cp_cfm_type) :: work
1361 :
1362 5953 : CALL timeset(routineN, handle)
1363 :
1364 5953 : IF (.NOT. PRESENT(workspace)) THEN
1365 5448 : CALL cp_cfm_create(work, matrix%matrix_struct)
1366 : ELSE
1367 505 : work = workspace
1368 : END IF
1369 :
1370 5953 : myuplo = 'U'
1371 5953 : IF (PRESENT(uplo)) myuplo = uplo
1372 :
1373 : ! get info of fm_mat_Q
1374 : CALL cp_cfm_get_info(matrix=matrix, &
1375 : nrow_local=nrow_local, &
1376 : ncol_local=ncol_local, &
1377 : row_indices=row_indices, &
1378 5953 : col_indices=col_indices)
1379 :
1380 247843 : DO jjB = 1, ncol_local
1381 241890 : j_global = col_indices(jjB)
1382 8307838 : DO iiB = 1, nrow_local
1383 8059995 : i_global = row_indices(iiB)
1384 8301885 : IF (MERGE(j_global < i_global, j_global > i_global, (myuplo == "U") .OR. (myuplo == "u"))) THEN
1385 3967260 : matrix%local_data(iiB, jjB) = z_zero
1386 4092735 : ELSE IF (j_global == i_global) THEN
1387 125475 : matrix%local_data(iiB, jjB) = matrix%local_data(iiB, jjB)/(2.0_dp, 0.0_dp)
1388 : END IF
1389 : END DO
1390 : END DO
1391 :
1392 5953 : CALL cp_cfm_transpose(matrix, 'C', work)
1393 :
1394 5953 : CALL cp_cfm_scale_and_add(z_one, matrix, z_one, work)
1395 :
1396 5953 : IF (.NOT. PRESENT(workspace)) THEN
1397 5448 : CALL cp_cfm_release(work)
1398 : END IF
1399 :
1400 5953 : CALL timestop(handle)
1401 :
1402 5953 : END SUBROUTINE cp_cfm_uplo_to_full
1403 :
1404 : ! **************************************************************************************************
1405 : !> \brief find the norm of each column norm_{j}= sqrt( \sum_{i} A_{ij}*conjg(A_{ij}) )
1406 : !> Complex-valued mirror of cp_fm_vectorsnorm.
1407 : !> \param matrix ...
1408 : !> \param norm_array ...
1409 : ! **************************************************************************************************
1410 42305 : SUBROUTINE cp_cfm_vectorsnorm(matrix, norm_array)
1411 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1412 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: norm_array
1413 :
1414 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_cfm_vectorsnorm'
1415 :
1416 : INTEGER :: handle, i, j, ncol_local, nrow_local
1417 42305 : INTEGER, DIMENSION(:), POINTER :: col_indices
1418 :
1419 42305 : CALL timeset(routineN, handle)
1420 :
1421 : CALL cp_cfm_get_info(matrix, col_indices=col_indices, nrow_local=nrow_local, &
1422 42305 : ncol_local=ncol_local)
1423 :
1424 : ! the efficiency could be improved by making use of the row-col distribution of scalapack
1425 696190 : norm_array = 0.0_dp
1426 696190 : DO j = 1, ncol_local
1427 20920002 : DO i = 1, nrow_local
1428 : norm_array(col_indices(j)) = norm_array(col_indices(j)) + &
1429 : REAL(matrix%local_data(i, j), KIND=dp)**2 + &
1430 20877697 : AIMAG(matrix%local_data(i, j))**2
1431 : END DO
1432 : END DO
1433 1350075 : CALL matrix%matrix_struct%para_env%sum(norm_array)
1434 696190 : norm_array = SQRT(norm_array)
1435 :
1436 42305 : CALL timestop(handle)
1437 42305 : END SUBROUTINE cp_cfm_vectorsnorm
1438 :
1439 : ! **************************************************************************************************
1440 : !> \brief returns the diagonal of a complex full matrix: diag(i)= A_{ii}.
1441 : !> Each diagonal entry is owned by one process. The sum over the
1442 : !> process grid collects the entries.
1443 : !> \param matrix ...
1444 : !> \param diag ...
1445 : ! **************************************************************************************************
1446 14955 : SUBROUTINE cp_cfm_get_diag(matrix, diag)
1447 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1448 : COMPLEX(KIND=dp), DIMENSION(:), INTENT(OUT) :: diag
1449 :
1450 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_cfm_get_diag'
1451 :
1452 : INTEGER :: handle, i, j, ncol_local, nrow_local
1453 14955 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1454 :
1455 14955 : CALL timeset(routineN, handle)
1456 :
1457 : CALL cp_cfm_get_info(matrix, col_indices=col_indices, row_indices=row_indices, &
1458 14955 : nrow_local=nrow_local, ncol_local=ncol_local)
1459 :
1460 235926 : diag = z_zero
1461 235926 : DO j = 1, ncol_local
1462 2762746 : DO i = 1, nrow_local
1463 2747791 : IF (row_indices(i) == col_indices(j)) THEN
1464 133320 : diag(col_indices(j)) = matrix%local_data(i, j)
1465 : END IF
1466 : END DO
1467 : END DO
1468 456897 : CALL matrix%matrix_struct%para_env%sum(diag)
1469 :
1470 14955 : CALL timestop(handle)
1471 14955 : END SUBROUTINE cp_cfm_get_diag
1472 :
1473 : END MODULE cp_cfm_basic_linalg
|
