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