Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 8232 : 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 204 : 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 204 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b, c
159 : INTEGER :: handle, icol_local, irow_local, mypcol, &
160 : myprow, ncol_local, nrow_local
161 :
162 204 : CALL timeset(routineN, handle)
163 :
164 204 : myprow = matrix_a%matrix_struct%context%mepos(1)
165 204 : mypcol = matrix_a%matrix_struct%context%mepos(2)
166 :
167 204 : a => matrix_a%local_data
168 204 : b => matrix_b%local_data
169 204 : c => matrix_c%local_data
170 :
171 204 : nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
172 204 : ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
173 :
174 1020 : DO icol_local = 1, ncol_local
175 2652 : DO irow_local = 1, nrow_local
176 2448 : c(irow_local, icol_local) = a(irow_local, icol_local)*b(irow_local, icol_local)
177 : END DO
178 : END DO
179 :
180 204 : CALL timestop(handle)
181 :
182 204 : 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 236168 : 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 236168 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b
249 : INTEGER :: handle, icol_local, irow_local, mypcol, &
250 : myprow, ncol_local, nrow_local
251 :
252 236168 : CALL timeset(routineN, handle)
253 :
254 236168 : my_beta = z_zero
255 236168 : IF (PRESENT(beta)) my_beta = beta
256 236168 : NULLIFY (a, b)
257 :
258 : ! to do: use dscal,dcopy,daxp
259 236168 : myprow = matrix_a%matrix_struct%context%mepos(1)
260 236168 : mypcol = matrix_a%matrix_struct%context%mepos(2)
261 :
262 236168 : nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
263 236168 : ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
264 :
265 236168 : a => matrix_a%local_data
266 :
267 236168 : IF (my_beta == z_zero) THEN
268 :
269 9064 : IF (alpha == z_zero) THEN
270 0 : a(:, :) = z_zero
271 9064 : ELSE IF (alpha == z_one) THEN
272 9064 : CALL timestop(handle)
273 9064 : RETURN
274 : ELSE
275 0 : a(:, :) = alpha*a(:, :)
276 : END IF
277 :
278 : ELSE
279 227104 : CPASSERT(PRESENT(matrix_b))
280 227104 : IF (matrix_a%matrix_struct%context /= matrix_b%matrix_struct%context) &
281 0 : CPABORT("matrixes must be in the same blacs context")
282 :
283 227104 : IF (cp_fm_struct_equivalent(matrix_a%matrix_struct, &
284 : matrix_b%matrix_struct)) THEN
285 :
286 227104 : b => matrix_b%local_data
287 :
288 227104 : 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 227104 : ELSE IF (alpha == z_one) THEN
305 223944 : IF (my_beta == z_one) THEN
306 : !a(:, :) = a(:, :)+b(:, :)
307 1903721 : DO icol_local = 1, ncol_local
308 27498250 : DO irow_local = 1, nrow_local
309 27334509 : 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 778273 : DO icol_local = 1, ncol_local
315 11241766 : DO irow_local = 1, nrow_local
316 11181563 : 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 38304 : DO icol_local = 1, ncol_local
323 782648 : DO irow_local = 1, nrow_local
324 779488 : 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 : #if defined(__parallel)
330 0 : CPABORT("to do (pdscal,pdcopy,pdaxpy)")
331 : #else
332 : CPABORT("")
333 : #endif
334 : END IF
335 : END IF
336 227104 : CALL timestop(handle)
337 236168 : 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 136866 : 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 136866 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
359 : INTEGER :: handle, icol_local, irow_local, mypcol, &
360 : myprow, ncol_local, nrow_local
361 136866 : REAL(kind=dp), DIMENSION(:, :), POINTER :: b
362 :
363 136866 : CALL timeset(routineN, handle)
364 :
365 136866 : NULLIFY (a, b)
366 :
367 136866 : myprow = matrix_a%matrix_struct%context%mepos(1)
368 136866 : mypcol = matrix_a%matrix_struct%context%mepos(2)
369 :
370 136866 : nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
371 136866 : ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
372 :
373 136866 : a => matrix_a%local_data
374 :
375 136866 : 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 136866 : IF (matrix_a%matrix_struct%context /= matrix_b%matrix_struct%context) &
388 0 : CPABORT("matrices must be in the same blacs context")
389 :
390 136866 : IF (cp_fm_struct_equivalent(matrix_a%matrix_struct, &
391 : matrix_b%matrix_struct)) THEN
392 :
393 136866 : b => matrix_b%local_data
394 :
395 136866 : IF (alpha == z_zero) THEN
396 50510 : IF (beta == z_one) THEN
397 : !a(:, :) = b(:, :)
398 1717002 : DO icol_local = 1, ncol_local
399 82207992 : DO irow_local = 1, nrow_local
400 82157506 : a(irow_local, icol_local) = b(irow_local, icol_local)
401 : END DO
402 : END DO
403 : ELSE
404 : !a(:, :) = beta*b(:, :)
405 456 : DO icol_local = 1, ncol_local
406 4344 : DO irow_local = 1, nrow_local
407 4320 : a(irow_local, icol_local) = beta*b(irow_local, icol_local)
408 : END DO
409 : END DO
410 : END IF
411 86356 : ELSE IF (alpha == z_one) THEN
412 55835 : IF (beta == z_one) THEN
413 : !a(:, :) = a(:, :)+b(:, :)
414 80449 : DO icol_local = 1, ncol_local
415 1438473 : DO irow_local = 1, nrow_local
416 1435080 : 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 1760566 : DO icol_local = 1, ncol_local
422 82715628 : DO irow_local = 1, nrow_local
423 82663186 : 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 346801 : DO icol_local = 1, ncol_local
430 5761521 : DO irow_local = 1, nrow_local
431 5731000 : 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 : #if defined(__parallel)
437 0 : CPABORT("to do (pdscal,pdcopy,pdaxpy)")
438 : #else
439 : CPABORT("")
440 : #endif
441 : END IF
442 : END IF
443 136866 : CALL timestop(handle)
444 136866 : END SUBROUTINE cp_cfm_scale_and_add_fm
445 :
446 : ! **************************************************************************************************
447 : !> \brief Computes LU decomposition of a given matrix.
448 : !> \param matrix_a full matrix
449 : !> \param determinant determinant
450 : !> \date 11.06.2001
451 : !> \author Matthias Krack
452 : !> \version 1.0
453 : !> \note
454 : !> The actual purpose right now is to efficiently compute the determinant of a given matrix.
455 : !> The original content of the matrix is destroyed.
456 : ! **************************************************************************************************
457 0 : SUBROUTINE cp_cfm_lu_decompose(matrix_a, determinant)
458 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
459 : COMPLEX(kind=dp), INTENT(out) :: determinant
460 :
461 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_lu_decompose'
462 :
463 0 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
464 : INTEGER :: counter, handle, info, irow, nrow_global
465 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
466 :
467 : #if defined(__parallel)
468 : INTEGER :: icol, ncol_local, nrow_local
469 : INTEGER, DIMENSION(9) :: desca
470 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
471 : #else
472 : INTEGER :: lda
473 : #endif
474 :
475 0 : CALL timeset(routineN, handle)
476 :
477 0 : nrow_global = matrix_a%matrix_struct%nrow_global
478 0 : a => matrix_a%local_data
479 :
480 0 : ALLOCATE (ipivot(nrow_global))
481 : #if defined(__parallel)
482 : CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, &
483 0 : row_indices=row_indices, col_indices=col_indices)
484 :
485 0 : desca(:) = matrix_a%matrix_struct%descriptor(:)
486 0 : CALL pzgetrf(nrow_global, nrow_global, a(1, 1), 1, 1, desca, ipivot, info)
487 :
488 0 : counter = 0
489 0 : DO irow = 1, nrow_local
490 0 : IF (ipivot(irow) .NE. row_indices(irow)) counter = counter + 1
491 : END DO
492 :
493 0 : IF (MOD(counter, 2) == 0) THEN
494 0 : determinant = z_one
495 : ELSE
496 0 : determinant = -z_one
497 : END IF
498 :
499 : ! compute product of diagonal elements
500 : irow = 1
501 : icol = 1
502 0 : DO WHILE (irow <= nrow_local .AND. icol <= ncol_local)
503 0 : IF (row_indices(irow) < col_indices(icol)) THEN
504 0 : irow = irow + 1
505 0 : ELSE IF (row_indices(irow) > col_indices(icol)) THEN
506 0 : icol = icol + 1
507 : ELSE ! diagonal element
508 0 : determinant = determinant*a(irow, icol)
509 0 : irow = irow + 1
510 0 : icol = icol + 1
511 : END IF
512 : END DO
513 0 : CALL matrix_a%matrix_struct%para_env%prod(determinant)
514 : #else
515 : lda = SIZE(a, 1)
516 : CALL zgetrf(nrow_global, nrow_global, a(1, 1), lda, ipivot, info)
517 : counter = 0
518 : determinant = z_one
519 : DO irow = 1, nrow_global
520 : IF (ipivot(irow) .NE. irow) counter = counter + 1
521 : determinant = determinant*a(irow, irow)
522 : END DO
523 : IF (MOD(counter, 2) == 1) determinant = -1.0_dp*determinant
524 : #endif
525 :
526 : ! info is allowed to be zero
527 : ! this does just signal a zero diagonal element
528 0 : DEALLOCATE (ipivot)
529 :
530 0 : CALL timestop(handle)
531 0 : END SUBROUTINE cp_cfm_lu_decompose
532 :
533 : ! **************************************************************************************************
534 : !> \brief Performs one of the matrix-matrix operations:
535 : !> matrix_c = alpha * op1( matrix_a ) * op2( matrix_b ) + beta*matrix_c.
536 : !> \param transa form of op1( matrix_a ):
537 : !> op1( matrix_a ) = matrix_a, when transa == 'N' ,
538 : !> op1( matrix_a ) = matrix_a^T, when transa == 'T' ,
539 : !> op1( matrix_a ) = matrix_a^H, when transa == 'C' ,
540 : !> \param transb form of op2( matrix_b )
541 : !> \param m number of rows of the matrix op1( matrix_a )
542 : !> \param n number of columns of the matrix op2( matrix_b )
543 : !> \param k number of columns of the matrix op1( matrix_a ) as well as
544 : !> number of rows of the matrix op2( matrix_b )
545 : !> \param alpha scale factor
546 : !> \param matrix_a matrix A
547 : !> \param matrix_b matrix B
548 : !> \param beta scale factor
549 : !> \param matrix_c matrix C
550 : !> \param a_first_col (optional) the first column of the matrix_a to multiply
551 : !> \param a_first_row (optional) the first row of the matrix_a to multiply
552 : !> \param b_first_col (optional) the first column of the matrix_b to multiply
553 : !> \param b_first_row (optional) the first row of the matrix_b to multiply
554 : !> \param c_first_col (optional) the first column of the matrix_c
555 : !> \param c_first_row (optional) the first row of the matrix_c
556 : !> \date 07.06.2001
557 : !> \author Matthias Krack
558 : !> \version 1.0
559 : ! **************************************************************************************************
560 76128 : SUBROUTINE cp_cfm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, &
561 : matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, &
562 : c_first_row)
563 : CHARACTER(len=1), INTENT(in) :: transa, transb
564 : INTEGER, INTENT(in) :: m, n, k
565 : COMPLEX(kind=dp), INTENT(in) :: alpha
566 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, matrix_b
567 : COMPLEX(kind=dp), INTENT(in) :: beta
568 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_c
569 : INTEGER, INTENT(in), OPTIONAL :: a_first_col, a_first_row, b_first_col, &
570 : b_first_row, c_first_col, c_first_row
571 :
572 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_gemm'
573 :
574 76128 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b, c
575 : INTEGER :: handle, i_a, i_b, i_c, j_a, j_b, j_c
576 : #if defined(__parallel)
577 : INTEGER, DIMENSION(9) :: desca, descb, descc
578 : #else
579 : INTEGER :: lda, ldb, ldc
580 : #endif
581 :
582 76128 : CALL timeset(routineN, handle)
583 76128 : a => matrix_a%local_data
584 76128 : b => matrix_b%local_data
585 76128 : c => matrix_c%local_data
586 :
587 76128 : i_a = 1
588 76128 : IF (PRESENT(a_first_row)) i_a = a_first_row
589 :
590 76128 : j_a = 1
591 76128 : IF (PRESENT(a_first_col)) j_a = a_first_col
592 :
593 76128 : i_b = 1
594 76128 : IF (PRESENT(b_first_row)) i_b = b_first_row
595 :
596 76128 : j_b = 1
597 76128 : IF (PRESENT(b_first_col)) j_b = b_first_col
598 :
599 76128 : i_c = 1
600 76128 : IF (PRESENT(c_first_row)) i_c = c_first_row
601 :
602 76128 : j_c = 1
603 76128 : IF (PRESENT(c_first_col)) j_c = c_first_col
604 :
605 : #if defined(__parallel)
606 761280 : desca(:) = matrix_a%matrix_struct%descriptor(:)
607 761280 : descb(:) = matrix_b%matrix_struct%descriptor(:)
608 761280 : descc(:) = matrix_c%matrix_struct%descriptor(:)
609 :
610 : CALL pzgemm(transa, transb, m, n, k, alpha, a(1, 1), i_a, j_a, desca, &
611 76128 : b(1, 1), i_b, j_b, descb, beta, c(1, 1), i_c, j_c, descc)
612 : #else
613 : lda = SIZE(a, 1)
614 : ldb = SIZE(b, 1)
615 : ldc = SIZE(c, 1)
616 :
617 : ! consider zgemm3m
618 : CALL zgemm(transa, transb, m, n, k, alpha, a(i_a, j_a), &
619 : lda, b(i_b, j_b), ldb, beta, c(i_c, j_c), ldc)
620 : #endif
621 76128 : CALL timestop(handle)
622 76128 : END SUBROUTINE cp_cfm_gemm
623 :
624 : ! **************************************************************************************************
625 : !> \brief Scales columns of the full matrix by corresponding factors.
626 : !> \param matrix_a matrix to scale
627 : !> \param scaling scale factors for every column. The actual number of scaled columns is
628 : !> limited by the number of scale factors given or by the actual number of columns
629 : !> whichever is smaller.
630 : !> \author Joost VandeVondele
631 : ! **************************************************************************************************
632 5816 : SUBROUTINE cp_cfm_column_scale(matrix_a, scaling)
633 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
634 : COMPLEX(kind=dp), DIMENSION(:), INTENT(in) :: scaling
635 :
636 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_column_scale'
637 :
638 5816 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
639 : INTEGER :: handle, icol_local, ncol_local, &
640 : nrow_local
641 : #if defined(__parallel)
642 5816 : INTEGER, DIMENSION(:), POINTER :: col_indices
643 : #endif
644 :
645 5816 : CALL timeset(routineN, handle)
646 :
647 5816 : a => matrix_a%local_data
648 :
649 : #if defined(__parallel)
650 5816 : CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, col_indices=col_indices)
651 5816 : ncol_local = MIN(ncol_local, SIZE(scaling))
652 :
653 284404 : DO icol_local = 1, ncol_local
654 284404 : CALL zscal(nrow_local, scaling(col_indices(icol_local)), a(1, icol_local), 1)
655 : END DO
656 : #else
657 : nrow_local = SIZE(a, 1)
658 : ncol_local = MIN(SIZE(a, 2), SIZE(scaling))
659 :
660 : DO icol_local = 1, ncol_local
661 : CALL zscal(nrow_local, scaling(icol_local), a(1, icol_local), 1)
662 : END DO
663 : #endif
664 :
665 5816 : CALL timestop(handle)
666 5816 : END SUBROUTINE cp_cfm_column_scale
667 :
668 : ! **************************************************************************************************
669 : !> \brief Scales a complex matrix by a real number.
670 : !> matrix_a = alpha * matrix_b
671 : !> \param alpha scale factor
672 : !> \param matrix_a complex matrix to scale
673 : ! **************************************************************************************************
674 8518 : SUBROUTINE cp_cfm_dscale(alpha, matrix_a)
675 : REAL(kind=dp), INTENT(in) :: alpha
676 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
677 :
678 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_dscale'
679 :
680 8518 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
681 : INTEGER :: handle
682 :
683 8518 : CALL timeset(routineN, handle)
684 :
685 8518 : NULLIFY (a)
686 :
687 8518 : a => matrix_a%local_data
688 :
689 25554 : CALL zdscal(SIZE(a), alpha, a(1, 1), 1)
690 :
691 8518 : CALL timestop(handle)
692 8518 : END SUBROUTINE cp_cfm_dscale
693 :
694 : ! **************************************************************************************************
695 : !> \brief Scales a complex matrix by a complex number.
696 : !> matrix_a = alpha * matrix_b
697 : !> \param alpha scale factor
698 : !> \param matrix_a complex matrix to scale
699 : !> \note
700 : !> use cp_fm_set_all to zero (avoids problems with nan)
701 : ! **************************************************************************************************
702 16150 : SUBROUTINE cp_cfm_zscale(alpha, matrix_a)
703 : COMPLEX(kind=dp), INTENT(IN) :: alpha
704 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
705 :
706 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_zscale'
707 :
708 16150 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
709 : INTEGER :: handle, size_a
710 :
711 16150 : CALL timeset(routineN, handle)
712 :
713 16150 : NULLIFY (a)
714 :
715 16150 : a => matrix_a%local_data
716 16150 : size_a = SIZE(a, 1)*SIZE(a, 2)
717 :
718 16150 : CALL zscal(size_a, alpha, a(1, 1), 1)
719 :
720 16150 : CALL timestop(handle)
721 16150 : END SUBROUTINE cp_cfm_zscale
722 :
723 : ! **************************************************************************************************
724 : !> \brief Solve the system of linear equations A*b=A_general using LU decomposition.
725 : !> Pay attention that both matrices are overwritten on exit and that
726 : !> the result is stored into the matrix 'general_a'.
727 : !> \param matrix_a matrix A (overwritten on exit)
728 : !> \param general_a (input) matrix A_general, (output) matrix B
729 : !> \param determinant (optional) determinant
730 : !> \author Florian Schiffmann
731 : ! **************************************************************************************************
732 7644 : SUBROUTINE cp_cfm_solve(matrix_a, general_a, determinant)
733 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, general_a
734 : COMPLEX(kind=dp), OPTIONAL :: determinant
735 :
736 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_solve'
737 :
738 7644 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, a_general
739 : INTEGER :: counter, handle, info, irow, nrow_global
740 7644 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
741 :
742 : #if defined(__parallel)
743 : INTEGER :: icol, ncol_local, nrow_local
744 : INTEGER, DIMENSION(9) :: desca, descb
745 7644 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
746 : #else
747 : INTEGER :: lda, ldb
748 : #endif
749 :
750 7644 : CALL timeset(routineN, handle)
751 :
752 7644 : a => matrix_a%local_data
753 7644 : a_general => general_a%local_data
754 7644 : nrow_global = matrix_a%matrix_struct%nrow_global
755 22932 : ALLOCATE (ipivot(nrow_global))
756 :
757 : #if defined(__parallel)
758 76440 : desca(:) = matrix_a%matrix_struct%descriptor(:)
759 76440 : descb(:) = general_a%matrix_struct%descriptor(:)
760 7644 : CALL pzgetrf(nrow_global, nrow_global, a(1, 1), 1, 1, desca, ipivot, info)
761 7644 : IF (PRESENT(determinant)) THEN
762 : CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, &
763 6356 : row_indices=row_indices, col_indices=col_indices)
764 :
765 6356 : counter = 0
766 19116 : DO irow = 1, nrow_local
767 19116 : IF (ipivot(irow) .NE. row_indices(irow)) counter = counter + 1
768 : END DO
769 :
770 6356 : IF (MOD(counter, 2) == 0) THEN
771 6346 : determinant = z_one
772 : ELSE
773 10 : determinant = -z_one
774 : END IF
775 :
776 : ! compute product of diagonal elements
777 : irow = 1
778 : icol = 1
779 28662 : DO WHILE (irow <= nrow_local .AND. icol <= ncol_local)
780 28662 : IF (row_indices(irow) < col_indices(icol)) THEN
781 0 : irow = irow + 1
782 22306 : ELSE IF (row_indices(irow) > col_indices(icol)) THEN
783 9546 : icol = icol + 1
784 : ELSE ! diagonal element
785 12760 : determinant = determinant*a(irow, icol)
786 12760 : irow = irow + 1
787 12760 : icol = icol + 1
788 : END IF
789 : END DO
790 6356 : CALL matrix_a%matrix_struct%para_env%prod(determinant)
791 : END IF
792 :
793 : CALL pzgetrs("N", nrow_global, nrow_global, a(1, 1), 1, 1, desca, &
794 7644 : ipivot, a_general(1, 1), 1, 1, descb, info)
795 : #else
796 : lda = SIZE(a, 1)
797 : ldb = SIZE(a_general, 1)
798 : CALL zgetrf(nrow_global, nrow_global, a(1, 1), lda, ipivot, info)
799 : IF (PRESENT(determinant)) THEN
800 : counter = 0
801 : determinant = z_one
802 : DO irow = 1, nrow_global
803 : IF (ipivot(irow) .NE. irow) counter = counter + 1
804 : determinant = determinant*a(irow, irow)
805 : END DO
806 : IF (MOD(counter, 2) == 1) determinant = -1.0_dp*determinant
807 : END IF
808 : CALL zgetrs("N", nrow_global, nrow_global, a(1, 1), lda, ipivot, a_general(1, 1), ldb, info)
809 : #endif
810 :
811 : ! info is allowed to be zero
812 : ! this does just signal a zero diagonal element
813 7644 : DEALLOCATE (ipivot)
814 7644 : CALL timestop(handle)
815 :
816 7644 : END SUBROUTINE cp_cfm_solve
817 :
818 : ! **************************************************************************************************
819 : !> \brief Inverts a matrix using LU decomposition. The input matrix will be overwritten.
820 : !> \param matrix input a general square non-singular matrix, outputs its inverse
821 : !> \param info_out optional, if present outputs the info from (p)zgetri
822 : !> \author Lianheng Tong
823 : ! **************************************************************************************************
824 49793 : SUBROUTINE cp_cfm_lu_invert(matrix, info_out)
825 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
826 : INTEGER, INTENT(out), OPTIONAL :: info_out
827 :
828 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_lu_invert'
829 :
830 49793 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
831 : COMPLEX(kind=dp), DIMENSION(1) :: work1
832 49793 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: mat
833 : INTEGER :: handle, info, lwork, nrows_global
834 49793 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
835 :
836 : #if defined(__parallel)
837 : INTEGER :: liwork
838 49793 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
839 : INTEGER, DIMENSION(1) :: iwork1
840 : INTEGER, DIMENSION(9) :: desca
841 : #else
842 : INTEGER :: lda
843 : #endif
844 :
845 49793 : CALL timeset(routineN, handle)
846 :
847 49793 : mat => matrix%local_data
848 49793 : nrows_global = matrix%matrix_struct%nrow_global
849 49793 : CPASSERT(nrows_global .EQ. matrix%matrix_struct%ncol_global)
850 149379 : ALLOCATE (ipivot(nrows_global))
851 :
852 : ! do LU decomposition
853 : #if defined(__parallel)
854 497930 : desca = matrix%matrix_struct%descriptor
855 : CALL pzgetrf(nrows_global, nrows_global, &
856 49793 : mat(1, 1), 1, 1, desca, ipivot, info)
857 : #else
858 : lda = SIZE(mat, 1)
859 : CALL zgetrf(nrows_global, nrows_global, &
860 : mat(1, 1), lda, ipivot, info)
861 : #endif
862 49793 : IF (info /= 0) THEN
863 0 : CALL cp_abort(__LOCATION__, "LU decomposition has failed")
864 : END IF
865 :
866 : ! do inversion
867 : #if defined(__parallel)
868 : CALL pzgetri(nrows_global, mat(1, 1), 1, 1, desca, &
869 49793 : ipivot, work1, -1, iwork1, -1, info)
870 49793 : lwork = INT(work1(1))
871 49793 : liwork = INT(iwork1(1))
872 149379 : ALLOCATE (work(lwork))
873 149379 : ALLOCATE (iwork(liwork))
874 : CALL pzgetri(nrows_global, mat(1, 1), 1, 1, desca, &
875 49793 : ipivot, work, lwork, iwork, liwork, info)
876 49793 : DEALLOCATE (iwork)
877 : #else
878 : CALL zgetri(nrows_global, mat(1, 1), lda, ipivot, work1, -1, info)
879 : lwork = INT(work1(1))
880 : ALLOCATE (work(lwork))
881 : CALL zgetri(nrows_global, mat(1, 1), lda, ipivot, work, lwork, info)
882 : #endif
883 49793 : DEALLOCATE (work)
884 49793 : DEALLOCATE (ipivot)
885 :
886 49793 : IF (PRESENT(info_out)) THEN
887 0 : info_out = info
888 : ELSE
889 49793 : IF (info /= 0) &
890 0 : CALL cp_abort(__LOCATION__, "LU inversion has failed")
891 : END IF
892 :
893 49793 : CALL timestop(handle)
894 :
895 49793 : END SUBROUTINE cp_cfm_lu_invert
896 :
897 : ! **************************************************************************************************
898 : !> \brief Returns the trace of matrix_a^T matrix_b, i.e
899 : !> sum_{i,j}(matrix_a(i,j)*matrix_b(i,j)) .
900 : !> \param matrix_a a complex matrix
901 : !> \param matrix_b another complex matrix
902 : !> \param trace value of the trace operator
903 : !> \par History
904 : !> * 09.2017 created [Sergey Chulkov]
905 : !> \author Sergey Chulkov
906 : !> \note
907 : !> Based on the subroutine cp_fm_trace(). Note the transposition of matrix_a!
908 : ! **************************************************************************************************
909 31315 : SUBROUTINE cp_cfm_trace(matrix_a, matrix_b, trace)
910 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, matrix_b
911 : COMPLEX(kind=dp), INTENT(out) :: trace
912 :
913 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_trace'
914 :
915 : INTEGER :: handle, mypcol, myprow, ncol_local, &
916 : npcol, nprow, nrow_local
917 : TYPE(cp_blacs_env_type), POINTER :: context
918 : TYPE(mp_comm_type) :: group
919 :
920 31315 : CALL timeset(routineN, handle)
921 :
922 31315 : context => matrix_a%matrix_struct%context
923 31315 : myprow = context%mepos(1)
924 31315 : mypcol = context%mepos(2)
925 31315 : nprow = context%num_pe(1)
926 31315 : npcol = context%num_pe(2)
927 :
928 31315 : group = matrix_a%matrix_struct%para_env
929 :
930 31315 : nrow_local = MIN(matrix_a%matrix_struct%nrow_locals(myprow), matrix_b%matrix_struct%nrow_locals(myprow))
931 31315 : ncol_local = MIN(matrix_a%matrix_struct%ncol_locals(mypcol), matrix_b%matrix_struct%ncol_locals(mypcol))
932 :
933 : ! compute an accurate dot-product
934 : trace = accurate_dot_product(matrix_a%local_data(1:nrow_local, 1:ncol_local), &
935 31315 : matrix_b%local_data(1:nrow_local, 1:ncol_local))
936 :
937 31315 : CALL group%sum(trace)
938 :
939 31315 : CALL timestop(handle)
940 :
941 31315 : END SUBROUTINE cp_cfm_trace
942 :
943 : ! **************************************************************************************************
944 : !> \brief Multiplies in place by a triangular matrix:
945 : !> matrix_b = alpha op(triangular_matrix) matrix_b
946 : !> or (if side='R')
947 : !> matrix_b = alpha matrix_b op(triangular_matrix)
948 : !> op(triangular_matrix) is:
949 : !> triangular_matrix (if transa="N" and invert_tr=.false.)
950 : !> triangular_matrix^T (if transa="T" and invert_tr=.false.)
951 : !> triangular_matrix^H (if transa="C" and invert_tr=.false.)
952 : !> triangular_matrix^(-1) (if transa="N" and invert_tr=.true.)
953 : !> triangular_matrix^(-T) (if transa="T" and invert_tr=.true.)
954 : !> triangular_matrix^(-H) (if transa="C" and invert_tr=.true.)
955 : !> \param triangular_matrix the triangular matrix that multiplies the other
956 : !> \param matrix_b the matrix that gets multiplied and stores the result
957 : !> \param side on which side of matrix_b stays op(triangular_matrix)
958 : !> (defaults to 'L')
959 : !> \param transa_tr ...
960 : !> \param invert_tr if the triangular matrix should be inverted
961 : !> (defaults to false)
962 : !> \param uplo_tr if triangular_matrix is stored in the upper ('U') or
963 : !> lower ('L') triangle (defaults to 'U')
964 : !> \param unit_diag_tr if the diagonal elements of triangular_matrix should
965 : !> be assumed to be 1 (defaults to false)
966 : !> \param n_rows the number of rows of the result (defaults to
967 : !> size(matrix_b,1))
968 : !> \param n_cols the number of columns of the result (defaults to
969 : !> size(matrix_b,2))
970 : !> \param alpha ...
971 : !> \par History
972 : !> 08.2002 created [fawzi]
973 : !> \author Fawzi Mohamed
974 : !> \note
975 : !> needs an mpi env
976 : ! **************************************************************************************************
977 101328 : SUBROUTINE cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, &
978 : transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, &
979 : alpha)
980 : TYPE(cp_cfm_type), INTENT(IN) :: triangular_matrix, matrix_b
981 : CHARACTER, INTENT(in), OPTIONAL :: side, transa_tr
982 : LOGICAL, INTENT(in), OPTIONAL :: invert_tr
983 : CHARACTER, INTENT(in), OPTIONAL :: uplo_tr
984 : LOGICAL, INTENT(in), OPTIONAL :: unit_diag_tr
985 : INTEGER, INTENT(in), OPTIONAL :: n_rows, n_cols
986 : COMPLEX(kind=dp), INTENT(in), OPTIONAL :: alpha
987 :
988 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_triangular_multiply'
989 :
990 : CHARACTER :: side_char, transa, unit_diag, uplo
991 : COMPLEX(kind=dp) :: al
992 : INTEGER :: handle, m, n
993 : LOGICAL :: invert
994 :
995 50664 : CALL timeset(routineN, handle)
996 50664 : side_char = 'L'
997 50664 : unit_diag = 'N'
998 50664 : uplo = 'U'
999 50664 : transa = 'N'
1000 50664 : invert = .FALSE.
1001 50664 : al = CMPLX(1.0_dp, 0.0_dp, dp)
1002 50664 : CALL cp_cfm_get_info(matrix_b, nrow_global=m, ncol_global=n)
1003 50664 : IF (PRESENT(side)) side_char = side
1004 50664 : IF (PRESENT(invert_tr)) invert = invert_tr
1005 50664 : IF (PRESENT(uplo_tr)) uplo = uplo_tr
1006 50664 : IF (PRESENT(unit_diag_tr)) THEN
1007 0 : IF (unit_diag_tr) THEN
1008 0 : unit_diag = 'U'
1009 : ELSE
1010 : unit_diag = 'N'
1011 : END IF
1012 : END IF
1013 50664 : IF (PRESENT(transa_tr)) transa = transa_tr
1014 50664 : IF (PRESENT(alpha)) al = alpha
1015 50664 : IF (PRESENT(n_rows)) m = n_rows
1016 50664 : IF (PRESENT(n_cols)) n = n_cols
1017 :
1018 50664 : IF (invert) THEN
1019 :
1020 : #if defined(__parallel)
1021 : CALL pztrsm(side_char, uplo, transa, unit_diag, m, n, al, &
1022 : triangular_matrix%local_data(1, 1), 1, 1, &
1023 : triangular_matrix%matrix_struct%descriptor, &
1024 : matrix_b%local_data(1, 1), 1, 1, &
1025 534 : matrix_b%matrix_struct%descriptor(1))
1026 : #else
1027 : CALL ztrsm(side_char, uplo, transa, unit_diag, m, n, al, &
1028 : triangular_matrix%local_data(1, 1), &
1029 : SIZE(triangular_matrix%local_data, 1), &
1030 : matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1))
1031 : #endif
1032 :
1033 : ELSE
1034 :
1035 : #if defined(__parallel)
1036 : CALL pztrmm(side_char, uplo, transa, unit_diag, m, n, al, &
1037 : triangular_matrix%local_data(1, 1), 1, 1, &
1038 : triangular_matrix%matrix_struct%descriptor, &
1039 : matrix_b%local_data(1, 1), 1, 1, &
1040 50130 : matrix_b%matrix_struct%descriptor(1))
1041 : #else
1042 : CALL ztrmm(side_char, uplo, transa, unit_diag, m, n, al, &
1043 : triangular_matrix%local_data(1, 1), &
1044 : SIZE(triangular_matrix%local_data, 1), &
1045 : matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1))
1046 : #endif
1047 :
1048 : END IF
1049 :
1050 50664 : CALL timestop(handle)
1051 :
1052 50664 : END SUBROUTINE cp_cfm_triangular_multiply
1053 :
1054 : ! **************************************************************************************************
1055 : !> \brief Inverts a triangular matrix.
1056 : !> \param matrix_a ...
1057 : !> \param uplo ...
1058 : !> \param info_out ...
1059 : !> \author MI
1060 : ! **************************************************************************************************
1061 16710 : SUBROUTINE cp_cfm_triangular_invert(matrix_a, uplo, info_out)
1062 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
1063 : CHARACTER, INTENT(in), OPTIONAL :: uplo
1064 : INTEGER, INTENT(out), OPTIONAL :: info_out
1065 :
1066 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_triangular_invert'
1067 :
1068 : CHARACTER :: unit_diag, my_uplo
1069 : INTEGER :: handle, info, ncol_global
1070 : COMPLEX(kind=dp), DIMENSION(:, :), &
1071 16710 : POINTER :: a
1072 : #if defined(__parallel)
1073 : INTEGER, DIMENSION(9) :: desca
1074 : #endif
1075 :
1076 16710 : CALL timeset(routineN, handle)
1077 :
1078 16710 : unit_diag = 'N'
1079 16710 : my_uplo = 'U'
1080 16710 : IF (PRESENT(uplo)) my_uplo = uplo
1081 :
1082 16710 : ncol_global = matrix_a%matrix_struct%ncol_global
1083 :
1084 16710 : a => matrix_a%local_data
1085 :
1086 : #if defined(__parallel)
1087 167100 : desca(:) = matrix_a%matrix_struct%descriptor(:)
1088 16710 : CALL pztrtri(my_uplo, unit_diag, ncol_global, a(1, 1), 1, 1, desca, info)
1089 : #else
1090 : CALL ztrtri(my_uplo, unit_diag, ncol_global, a(1, 1), ncol_global, info)
1091 : #endif
1092 :
1093 16710 : IF (PRESENT(info_out)) THEN
1094 0 : info_out = info
1095 : ELSE
1096 16710 : IF (info /= 0) &
1097 : CALL cp_abort(__LOCATION__, &
1098 0 : "triangular invert failed: matrix is not positive definite or ill-conditioned")
1099 : END IF
1100 :
1101 16710 : CALL timestop(handle)
1102 16710 : END SUBROUTINE cp_cfm_triangular_invert
1103 :
1104 : ! **************************************************************************************************
1105 : !> \brief Transposes a BLACS distributed complex matrix.
1106 : !> \param matrix input matrix
1107 : !> \param trans 'T' for transpose, 'C' for Hermitian conjugate
1108 : !> \param matrixt output matrix
1109 : !> \author Lianheng Tong
1110 : ! **************************************************************************************************
1111 14506 : SUBROUTINE cp_cfm_transpose(matrix, trans, matrixt)
1112 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1113 : CHARACTER, INTENT(in) :: trans
1114 : TYPE(cp_cfm_type), INTENT(IN) :: matrixt
1115 :
1116 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_transpose'
1117 :
1118 14506 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa, cc
1119 : INTEGER :: handle, ncol_global, nrow_global
1120 : #if defined(__parallel)
1121 : INTEGER, DIMENSION(9) :: desca, descc
1122 : #elif !defined(__MKL)
1123 : INTEGER :: ii, jj
1124 : #endif
1125 :
1126 14506 : CALL timeset(routineN, handle)
1127 :
1128 14506 : nrow_global = matrix%matrix_struct%nrow_global
1129 14506 : ncol_global = matrix%matrix_struct%ncol_global
1130 :
1131 14506 : CPASSERT(matrixt%matrix_struct%nrow_global == ncol_global)
1132 14506 : CPASSERT(matrixt%matrix_struct%ncol_global == nrow_global)
1133 :
1134 14506 : aa => matrix%local_data
1135 14506 : cc => matrixt%local_data
1136 :
1137 : #if defined(__parallel)
1138 145060 : desca = matrix%matrix_struct%descriptor
1139 145060 : descc = matrixt%matrix_struct%descriptor
1140 6610 : SELECT CASE (trans)
1141 : CASE ('T')
1142 : CALL pztranu(nrow_global, ncol_global, &
1143 : z_one, aa(1, 1), 1, 1, desca, &
1144 6610 : z_zero, cc(1, 1), 1, 1, descc)
1145 : CASE ('C')
1146 : CALL pztranc(nrow_global, ncol_global, &
1147 : z_one, aa(1, 1), 1, 1, desca, &
1148 7896 : z_zero, cc(1, 1), 1, 1, descc)
1149 : CASE DEFAULT
1150 14506 : CPABORT("trans only accepts 'T' or 'C'")
1151 : END SELECT
1152 : #elif defined(__MKL)
1153 : CALL mkl_zomatcopy('C', trans, nrow_global, ncol_global, 1.0_dp, aa(1, 1), nrow_global, cc(1, 1), ncol_global)
1154 : #else
1155 : SELECT CASE (trans)
1156 : CASE ('T')
1157 : DO jj = 1, ncol_global
1158 : DO ii = 1, nrow_global
1159 : cc(ii, jj) = aa(jj, ii)
1160 : END DO
1161 : END DO
1162 : CASE ('C')
1163 : DO jj = 1, ncol_global
1164 : DO ii = 1, nrow_global
1165 : cc(ii, jj) = CONJG(aa(jj, ii))
1166 : END DO
1167 : END DO
1168 : CASE DEFAULT
1169 : CPABORT("trans only accepts 'T' or 'C'")
1170 : END SELECT
1171 : #endif
1172 :
1173 14506 : CALL timestop(handle)
1174 14506 : END SUBROUTINE cp_cfm_transpose
1175 :
1176 : ! **************************************************************************************************
1177 : !> \brief Norm of matrix using (p)zlange.
1178 : !> \param matrix input a general matrix
1179 : !> \param mode 'M' max abs element value,
1180 : !> '1' or 'O' one norm, i.e. maximum column sum,
1181 : !> 'I' infinity norm, i.e. maximum row sum,
1182 : !> 'F' or 'E' Frobenius norm, i.e. sqrt of sum of all squares of elements
1183 : !> \return the norm according to mode
1184 : !> \author Lianheng Tong
1185 : ! **************************************************************************************************
1186 104498 : FUNCTION cp_cfm_norm(matrix, mode) RESULT(res)
1187 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1188 : CHARACTER, INTENT(IN) :: mode
1189 : REAL(kind=dp) :: res
1190 :
1191 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_norm'
1192 :
1193 104498 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa
1194 : INTEGER :: handle, lwork, ncols, ncols_local, &
1195 : nrows, nrows_local
1196 104498 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
1197 :
1198 : #if defined(__parallel)
1199 : INTEGER, DIMENSION(9) :: desca
1200 : #else
1201 : INTEGER :: lda
1202 : #endif
1203 :
1204 104498 : CALL timeset(routineN, handle)
1205 :
1206 : CALL cp_cfm_get_info(matrix=matrix, &
1207 : nrow_global=nrows, &
1208 : ncol_global=ncols, &
1209 : nrow_local=nrows_local, &
1210 104498 : ncol_local=ncols_local)
1211 104498 : aa => matrix%local_data
1212 :
1213 : SELECT CASE (mode)
1214 : CASE ('M', 'm')
1215 0 : lwork = 1
1216 : CASE ('1', 'O', 'o')
1217 : #if defined(__parallel)
1218 0 : lwork = ncols_local
1219 : #else
1220 : lwork = 1
1221 : #endif
1222 : CASE ('I', 'i')
1223 : #if defined(__parallel)
1224 0 : lwork = nrows_local
1225 : #else
1226 : lwork = nrows
1227 : #endif
1228 : CASE ('F', 'f', 'E', 'e')
1229 0 : lwork = 1
1230 : CASE DEFAULT
1231 104498 : CPABORT("mode input is not valid")
1232 : END SELECT
1233 :
1234 313494 : ALLOCATE (work(lwork))
1235 :
1236 : #if defined(__parallel)
1237 1044980 : desca = matrix%matrix_struct%descriptor
1238 104498 : res = pzlange(mode, nrows, ncols, aa(1, 1), 1, 1, desca, work)
1239 : #else
1240 : lda = SIZE(aa, 1)
1241 : res = zlange(mode, nrows, ncols, aa(1, 1), lda, work)
1242 : #endif
1243 :
1244 104498 : DEALLOCATE (work)
1245 104498 : CALL timestop(handle)
1246 104498 : END FUNCTION cp_cfm_norm
1247 :
1248 : ! **************************************************************************************************
1249 : !> \brief Applies a planar rotation defined by cs and sn to the i'th and j'th rows.
1250 : !> \param matrix ...
1251 : !> \param irow ...
1252 : !> \param jrow ...
1253 : !> \param cs cosine of the rotation angle
1254 : !> \param sn sinus of the rotation angle
1255 : !> \author Ole Schuett
1256 : ! **************************************************************************************************
1257 0 : SUBROUTINE cp_cfm_rot_rows(matrix, irow, jrow, cs, sn)
1258 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1259 : INTEGER, INTENT(IN) :: irow, jrow
1260 : REAL(dp), INTENT(IN) :: cs, sn
1261 :
1262 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_rot_rows'
1263 : INTEGER :: handle, ncol
1264 : COMPLEX(KIND=dp) :: sn_cmplx
1265 :
1266 : #if defined(__parallel)
1267 : INTEGER :: info, lwork
1268 : INTEGER, DIMENSION(9) :: desc
1269 0 : REAL(dp), DIMENSION(:), ALLOCATABLE :: work
1270 : #endif
1271 0 : CALL timeset(routineN, handle)
1272 0 : CALL cp_cfm_get_info(matrix, ncol_global=ncol)
1273 0 : sn_cmplx = CMPLX(sn, 0.0_dp, dp)
1274 : #if defined(__parallel)
1275 0 : IF (1 /= matrix%matrix_struct%context%n_pid) THEN
1276 0 : lwork = 2*ncol + 1
1277 0 : ALLOCATE (work(lwork))
1278 0 : desc(:) = matrix%matrix_struct%descriptor(:)
1279 : CALL pzrot(ncol, &
1280 : matrix%local_data(1, 1), irow, 1, desc, ncol, &
1281 : matrix%local_data(1, 1), jrow, 1, desc, ncol, &
1282 0 : cs, sn_cmplx, work, lwork, info)
1283 0 : CPASSERT(info == 0)
1284 0 : DEALLOCATE (work)
1285 : ELSE
1286 : #endif
1287 0 : CALL zrot(ncol, matrix%local_data(irow, 1), ncol, matrix%local_data(jrow, 1), ncol, cs, sn_cmplx)
1288 : #if defined(__parallel)
1289 : END IF
1290 : #endif
1291 0 : CALL timestop(handle)
1292 0 : END SUBROUTINE cp_cfm_rot_rows
1293 :
1294 : ! **************************************************************************************************
1295 : !> \brief Applies a planar rotation defined by cs and sn to the i'th and j'th columnns.
1296 : !> \param matrix ...
1297 : !> \param icol ...
1298 : !> \param jcol ...
1299 : !> \param cs cosine of the rotation angle
1300 : !> \param sn sinus of the rotation angle
1301 : !> \author Ole Schuett
1302 : ! **************************************************************************************************
1303 0 : SUBROUTINE cp_cfm_rot_cols(matrix, icol, jcol, cs, sn)
1304 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1305 : INTEGER, INTENT(IN) :: icol, jcol
1306 : REAL(dp), INTENT(IN) :: cs, sn
1307 :
1308 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_rot_cols'
1309 : INTEGER :: handle, nrow
1310 : COMPLEX(KIND=dp) :: sn_cmplx
1311 :
1312 : #if defined(__parallel)
1313 : INTEGER :: info, lwork
1314 : INTEGER, DIMENSION(9) :: desc
1315 0 : REAL(dp), DIMENSION(:), ALLOCATABLE :: work
1316 : #endif
1317 0 : CALL timeset(routineN, handle)
1318 0 : CALL cp_cfm_get_info(matrix, nrow_global=nrow)
1319 0 : sn_cmplx = CMPLX(sn, 0.0_dp, dp)
1320 : #if defined(__parallel)
1321 0 : IF (1 /= matrix%matrix_struct%context%n_pid) THEN
1322 0 : lwork = 2*nrow + 1
1323 0 : ALLOCATE (work(lwork))
1324 0 : desc(:) = matrix%matrix_struct%descriptor(:)
1325 : CALL pzrot(nrow, &
1326 : matrix%local_data(1, 1), 1, icol, desc, 1, &
1327 : matrix%local_data(1, 1), 1, jcol, desc, 1, &
1328 0 : cs, sn_cmplx, work, lwork, info)
1329 0 : CPASSERT(info == 0)
1330 0 : DEALLOCATE (work)
1331 : ELSE
1332 : #endif
1333 0 : CALL zrot(nrow, matrix%local_data(1, icol), 1, matrix%local_data(1, jcol), 1, cs, sn_cmplx)
1334 : #if defined(__parallel)
1335 : END IF
1336 : #endif
1337 0 : CALL timestop(handle)
1338 0 : END SUBROUTINE cp_cfm_rot_cols
1339 :
1340 : ! **************************************************************************************************
1341 : !> \brief ...
1342 : !> \param matrix ...
1343 : !> \param workspace ...
1344 : !> \param uplo triangular format; defaults to 'U'
1345 : !> \par History
1346 : !> 12.2024 Added optional workspace as input [Rocco Meli]
1347 : !> \author Jan Wilhelm
1348 : ! **************************************************************************************************
1349 9636 : SUBROUTINE cp_cfm_uplo_to_full(matrix, workspace, uplo)
1350 :
1351 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1352 : TYPE(cp_cfm_type), INTENT(IN), OPTIONAL :: workspace
1353 : CHARACTER, INTENT(IN), OPTIONAL :: uplo
1354 :
1355 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_cfm_uplo_to_full'
1356 :
1357 : CHARACTER :: myuplo
1358 : INTEGER :: handle, i_global, iiB, j_global, jjB, &
1359 : ncol_local, nrow_local
1360 4818 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1361 : TYPE(cp_cfm_type) :: work
1362 :
1363 4818 : CALL timeset(routineN, handle)
1364 :
1365 4818 : IF (.NOT. PRESENT(workspace)) THEN
1366 4212 : CALL cp_cfm_create(work, matrix%matrix_struct)
1367 : ELSE
1368 606 : work = workspace
1369 : END IF
1370 :
1371 4818 : myuplo = 'U'
1372 4818 : IF (PRESENT(uplo)) myuplo = uplo
1373 :
1374 : ! get info of fm_mat_Q
1375 : CALL cp_cfm_get_info(matrix=matrix, &
1376 : nrow_local=nrow_local, &
1377 : ncol_local=ncol_local, &
1378 : row_indices=row_indices, &
1379 4818 : col_indices=col_indices)
1380 :
1381 200784 : DO jjB = 1, ncol_local
1382 195966 : j_global = col_indices(jjB)
1383 7187880 : DO iiB = 1, nrow_local
1384 6987096 : i_global = row_indices(iiB)
1385 7183062 : IF (MERGE(j_global < i_global, j_global > i_global, (myuplo == "U") .OR. (myuplo == "u"))) THEN
1386 3443190 : matrix%local_data(iiB, jjB) = z_zero
1387 3543906 : ELSE IF (j_global == i_global) THEN
1388 100716 : matrix%local_data(iiB, jjB) = matrix%local_data(iiB, jjB)/(2.0_dp, 0.0_dp)
1389 : END IF
1390 : END DO
1391 : END DO
1392 :
1393 4818 : CALL cp_cfm_transpose(matrix, 'C', work)
1394 :
1395 4818 : CALL cp_cfm_scale_and_add(z_one, matrix, z_one, work)
1396 :
1397 4818 : IF (.NOT. PRESENT(workspace)) THEN
1398 4212 : CALL cp_cfm_release(work)
1399 : END IF
1400 :
1401 4818 : CALL timestop(handle)
1402 :
1403 4818 : END SUBROUTINE cp_cfm_uplo_to_full
1404 :
1405 : END MODULE cp_cfm_basic_linalg
|
