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 OMP_LIB, ONLY: omp_get_max_threads
18 : USE cp_blacs_env, ONLY: cp_blacs_env_type
19 : USE cp_cfm_types, ONLY: cp_cfm_create,&
20 : cp_cfm_get_info,&
21 : cp_cfm_release,&
22 : cp_cfm_to_cfm,&
23 : cp_cfm_type
24 : USE cp_fm_struct, ONLY: cp_fm_struct_equivalent
25 : USE cp_fm_types, ONLY: cp_fm_type
26 : USE cp_log_handling, ONLY: cp_to_string
27 : USE kahan_sum, ONLY: accurate_dot_product
28 : USE kinds, ONLY: dp
29 : USE mathconstants, ONLY: z_one,&
30 : z_zero
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_add_on_diag, &
40 : cp_cfm_column_scale, &
41 : cp_cfm_gemm, &
42 : cp_cfm_get_diag, &
43 : cp_cfm_lu_decompose, &
44 : cp_cfm_lu_invert, &
45 : cp_cfm_norm, &
46 : cp_cfm_scale, &
47 : cp_cfm_scale_and_add, &
48 : cp_cfm_scale_and_add_fm, &
49 : cp_cfm_schur_product, &
50 : cp_cfm_schur_product_cc, &
51 : cp_cfm_solve, &
52 : cp_cfm_trace, &
53 : cp_cfm_transpose, &
54 : cp_cfm_triangular_invert, &
55 : cp_cfm_triangular_multiply, &
56 : cp_cfm_rot_rows, &
57 : cp_cfm_rot_cols, &
58 : cp_cfm_det, & ! determinant of a complex matrix with correct sign
59 : cp_cfm_uplo_to_full, &
60 : cp_cfm_vectorsnorm
61 :
62 : REAL(kind=dp), EXTERNAL :: zlange, pzlange
63 :
64 : INTERFACE cp_cfm_scale
65 : MODULE PROCEDURE cp_cfm_dscale, cp_cfm_zscale
66 : END INTERFACE cp_cfm_scale
67 :
68 : ! **************************************************************************************************
69 :
70 : CONTAINS
71 :
72 : ! **************************************************************************************************
73 : !> \brief Computes the determinant (with a correct sign even in parallel environment!) of a complex square matrix
74 : !> \param matrix_a ...
75 : !> \param det_a ...
76 : !> \author A. Sinyavskiy (andrey.sinyavskiy@chem.uzh.ch)
77 : ! **************************************************************************************************
78 1536 : SUBROUTINE cp_cfm_det(matrix_a, det_a)
79 :
80 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
81 : COMPLEX(KIND=dp), INTENT(OUT) :: det_a
82 : COMPLEX(KIND=dp) :: determinant
83 : TYPE(cp_cfm_type) :: matrix_lu
84 1536 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a
85 : INTEGER :: n, i, info, P
86 1536 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
87 1536 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: diag
88 :
89 : #if defined(__parallel)
90 : INTEGER :: myprow, nprow, npcol, nrow_local, irow_local, &
91 : mypcol, ncol_local, icol_local, j
92 : INTEGER, DIMENSION(9) :: desca
93 : #endif
94 :
95 : CALL cp_cfm_create(matrix=matrix_lu, &
96 : matrix_struct=matrix_a%matrix_struct, &
97 1536 : name="A_lu"//TRIM(ADJUSTL(cp_to_string(1)))//"MATRIX")
98 1536 : CALL cp_cfm_to_cfm(matrix_a, matrix_lu)
99 :
100 1536 : a => matrix_lu%local_data
101 1536 : n = matrix_lu%matrix_struct%nrow_global
102 4608 : ALLOCATE (ipivot(n))
103 1536 : ipivot(:) = 0
104 1536 : P = 0
105 4608 : ALLOCATE (diag(n))
106 8622 : diag(:) = 0.0_dp
107 : #if defined(__parallel)
108 : ! Use LU decomposition
109 15360 : desca(:) = matrix_lu%matrix_struct%descriptor(:)
110 1536 : CALL pzgetrf(n, n, a(1, 1), 1, 1, desca, ipivot, info)
111 1536 : myprow = matrix_lu%matrix_struct%context%mepos(1)
112 1536 : mypcol = matrix_lu%matrix_struct%context%mepos(2)
113 1536 : nprow = matrix_lu%matrix_struct%context%num_pe(1)
114 1536 : npcol = matrix_lu%matrix_struct%context%num_pe(2)
115 1536 : nrow_local = matrix_lu%matrix_struct%nrow_locals(myprow)
116 1536 : ncol_local = matrix_lu%matrix_struct%ncol_locals(mypcol)
117 :
118 5199 : DO irow_local = 1, nrow_local
119 3663 : i = matrix_lu%matrix_struct%row_indices(irow_local)
120 41874 : DO icol_local = 1, ncol_local
121 36675 : j = matrix_lu%matrix_struct%col_indices(icol_local)
122 40338 : IF (i == j) diag(i) = matrix_lu%local_data(irow_local, icol_local)
123 : END DO
124 : END DO
125 15708 : CALL matrix_lu%matrix_struct%para_env%sum(diag)
126 8622 : determinant = PRODUCT(diag)
127 5199 : DO irow_local = 1, nrow_local
128 3663 : i = matrix_lu%matrix_struct%row_indices(irow_local)
129 5199 : IF (ipivot(irow_local) /= i) P = P + 1
130 : END DO
131 1536 : CALL matrix_lu%matrix_struct%para_env%sum(P)
132 : ! very important fix
133 1536 : P = P/npcol
134 : #else
135 : CALL zgetrf(n, n, a(1, 1), n, ipivot, info)
136 : DO i = 1, n
137 : diag(i) = matrix_lu%local_data(i, i)
138 : END DO
139 : determinant = PRODUCT(diag)
140 : DO i = 1, n
141 : IF (ipivot(i) /= i) P = P + 1
142 : END DO
143 : #endif
144 1536 : DEALLOCATE (ipivot)
145 1536 : DEALLOCATE (diag)
146 1536 : CALL cp_cfm_release(matrix_lu)
147 1536 : det_a = determinant*(-2*MOD(P, 2) + 1.0_dp)
148 1536 : END SUBROUTINE cp_cfm_det
149 :
150 : ! **************************************************************************************************
151 : !> \brief Computes the element-wise (Schur) product of two matrices: C = A \circ B .
152 : !> \param matrix_a the first input matrix
153 : !> \param matrix_b the second input matrix
154 : !> \param matrix_c matrix to store the result
155 : ! **************************************************************************************************
156 154 : SUBROUTINE cp_cfm_schur_product(matrix_a, matrix_b, matrix_c)
157 :
158 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, matrix_b, matrix_c
159 :
160 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_schur_product'
161 :
162 154 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b, c
163 : INTEGER :: handle, icol_local, irow_local, mypcol, &
164 : myprow, ncol_local, nrow_local
165 :
166 154 : CALL timeset(routineN, handle)
167 :
168 154 : myprow = matrix_a%matrix_struct%context%mepos(1)
169 154 : mypcol = matrix_a%matrix_struct%context%mepos(2)
170 :
171 154 : a => matrix_a%local_data
172 154 : b => matrix_b%local_data
173 154 : c => matrix_c%local_data
174 :
175 154 : nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
176 154 : ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
177 :
178 462 : DO icol_local = 1, ncol_local
179 770 : DO irow_local = 1, nrow_local
180 616 : c(irow_local, icol_local) = a(irow_local, icol_local)*b(irow_local, icol_local)
181 : END DO
182 : END DO
183 :
184 154 : CALL timestop(handle)
185 :
186 154 : END SUBROUTINE cp_cfm_schur_product
187 :
188 : ! **************************************************************************************************
189 : !> \brief Computes the element-wise (Schur) product of two matrices: C = A \circ conjg(B) .
190 : !> \param matrix_a the first input matrix
191 : !> \param matrix_b the second input matrix
192 : !> \param matrix_c matrix to store the result
193 : ! **************************************************************************************************
194 760 : SUBROUTINE cp_cfm_schur_product_cc(matrix_a, matrix_b, matrix_c)
195 :
196 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, matrix_b, matrix_c
197 :
198 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_schur_product_cc'
199 :
200 760 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b, c
201 : INTEGER :: handle, icol_local, irow_local, mypcol, &
202 : myprow, ncol_local, nrow_local
203 :
204 760 : CALL timeset(routineN, handle)
205 :
206 760 : myprow = matrix_a%matrix_struct%context%mepos(1)
207 760 : mypcol = matrix_a%matrix_struct%context%mepos(2)
208 :
209 760 : a => matrix_a%local_data
210 760 : b => matrix_b%local_data
211 760 : c => matrix_c%local_data
212 :
213 760 : nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
214 760 : ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
215 :
216 5233 : DO icol_local = 1, ncol_local
217 190760 : DO irow_local = 1, nrow_local
218 190000 : c(irow_local, icol_local) = a(irow_local, icol_local)*CONJG(b(irow_local, icol_local))
219 : END DO
220 : END DO
221 :
222 760 : CALL timestop(handle)
223 :
224 760 : END SUBROUTINE cp_cfm_schur_product_cc
225 :
226 : ! **************************************************************************************************
227 : !> \brief Scale and add two BLACS matrices (a = alpha*a + beta*b).
228 : !> \param alpha ...
229 : !> \param matrix_a ...
230 : !> \param beta ...
231 : !> \param matrix_b ...
232 : !> \date 11.06.2001
233 : !> \author Matthias Krack
234 : !> \version 1.0
235 : !> \note
236 : !> Use explicit loops to avoid temporary arrays, as a compiler reasonably assumes that arrays
237 : !> matrix_a%local_data and matrix_b%local_data may overlap (they are referenced by pointers).
238 : !> In general case (alpha*a + beta*b) explicit loops appears to be up to two times more efficient
239 : !> than equivalent LAPACK calls (zscale, zaxpy). This is because using LAPACK calls implies
240 : !> two passes through each array, so data need to be retrieved twice if arrays are large
241 : !> enough to not fit into the processor's cache.
242 : ! **************************************************************************************************
243 664574 : SUBROUTINE cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
244 : COMPLEX(kind=dp), INTENT(in) :: alpha
245 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
246 : COMPLEX(kind=dp), INTENT(in), OPTIONAL :: beta
247 : TYPE(cp_cfm_type), INTENT(IN), OPTIONAL :: matrix_b
248 :
249 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_scale_and_add'
250 :
251 : COMPLEX(kind=dp) :: my_beta
252 664574 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b
253 : INTEGER :: handle, icol_local, irow_local, mypcol, &
254 : myprow, ncol_local, nrow_local
255 :
256 664574 : CALL timeset(routineN, handle)
257 :
258 664574 : my_beta = z_zero
259 664574 : IF (PRESENT(beta)) my_beta = beta
260 664574 : NULLIFY (a, b)
261 :
262 : ! to do: use dscal,dcopy,daxp
263 664574 : myprow = matrix_a%matrix_struct%context%mepos(1)
264 664574 : mypcol = matrix_a%matrix_struct%context%mepos(2)
265 :
266 664574 : nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
267 664574 : ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
268 :
269 664574 : a => matrix_a%local_data
270 :
271 664574 : IF (my_beta == z_zero) THEN
272 :
273 65984 : IF (alpha == z_zero) THEN
274 0 : a(:, :) = z_zero
275 65984 : ELSE IF (alpha == z_one) THEN
276 65984 : CALL timestop(handle)
277 65984 : RETURN
278 : ELSE
279 0 : a(:, :) = alpha*a(:, :)
280 : END IF
281 :
282 : ELSE
283 598590 : CPASSERT(PRESENT(matrix_b))
284 598590 : IF (matrix_a%matrix_struct%context /= matrix_b%matrix_struct%context) &
285 0 : CPABORT("matrixes must be in the same blacs context")
286 :
287 598590 : IF (cp_fm_struct_equivalent(matrix_a%matrix_struct, &
288 : matrix_b%matrix_struct)) THEN
289 :
290 598590 : b => matrix_b%local_data
291 :
292 598590 : IF (alpha == z_zero) THEN
293 20 : IF (my_beta == z_one) THEN
294 : !a(:, :) = b(:, :)
295 724 : DO icol_local = 1, ncol_local
296 49172 : DO irow_local = 1, nrow_local
297 49152 : a(irow_local, icol_local) = b(irow_local, icol_local)
298 : END DO
299 : END DO
300 : ELSE
301 : !a(:, :) = my_beta*b(:, :)
302 0 : DO icol_local = 1, ncol_local
303 0 : DO irow_local = 1, nrow_local
304 0 : a(irow_local, icol_local) = my_beta*b(irow_local, icol_local)
305 : END DO
306 : END DO
307 : END IF
308 598570 : ELSE IF (alpha == z_one) THEN
309 578021 : IF (my_beta == z_one) THEN
310 : !a(:, :) = a(:, :)+b(:, :)
311 4221307 : DO icol_local = 1, ncol_local
312 60727771 : DO irow_local = 1, nrow_local
313 60309708 : a(irow_local, icol_local) = a(irow_local, icol_local) + b(irow_local, icol_local)
314 : END DO
315 : END DO
316 : ELSE
317 : !a(:, :) = a(:, :)+my_beta*b(:, :)
318 2484807 : DO icol_local = 1, ncol_local
319 72666645 : DO irow_local = 1, nrow_local
320 72506687 : a(irow_local, icol_local) = a(irow_local, icol_local) + my_beta*b(irow_local, icol_local)
321 : END DO
322 : END DO
323 : END IF
324 : ELSE
325 : !a(:, :) = alpha*a(:, :)+my_beta*b(:, :)
326 305196 : DO icol_local = 1, ncol_local
327 8082612 : DO irow_local = 1, nrow_local
328 8062063 : a(irow_local, icol_local) = alpha*a(irow_local, icol_local) + my_beta*b(irow_local, icol_local)
329 : END DO
330 : END DO
331 : END IF
332 : ELSE
333 : CALL cp_abort(__LOCATION__, &
334 : "cp_cfm_scale_and_add is not yet implemented for cases "// &
335 0 : "where input two matrix structures are not equivalent")
336 : END IF
337 : END IF
338 598590 : CALL timestop(handle)
339 664574 : END SUBROUTINE cp_cfm_scale_and_add
340 :
341 : ! **************************************************************************************************
342 : !> \brief Scale and add two BLACS matrices (a = alpha*a + beta*b).
343 : !> where b is a real matrix (adapted from cp_cfm_scale_and_add).
344 : !> \param alpha ...
345 : !> \param matrix_a ...
346 : !> \param beta ...
347 : !> \param matrix_b ...
348 : !> \date 01.08.2014
349 : !> \author JGH
350 : !> \version 1.0
351 : ! **************************************************************************************************
352 96016 : SUBROUTINE cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
353 : COMPLEX(kind=dp), INTENT(in) :: alpha
354 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
355 : COMPLEX(kind=dp), INTENT(in) :: beta
356 : TYPE(cp_fm_type), INTENT(IN) :: matrix_b
357 :
358 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_scale_and_add_fm'
359 :
360 96016 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
361 : INTEGER :: handle, icol_local, irow_local, mypcol, &
362 : myprow, ncol_local, nrow_local
363 96016 : REAL(kind=dp), DIMENSION(:, :), POINTER :: b
364 :
365 96016 : CALL timeset(routineN, handle)
366 :
367 96016 : NULLIFY (a, b)
368 :
369 96016 : myprow = matrix_a%matrix_struct%context%mepos(1)
370 96016 : mypcol = matrix_a%matrix_struct%context%mepos(2)
371 :
372 96016 : nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
373 96016 : ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
374 :
375 96016 : a => matrix_a%local_data
376 :
377 96016 : IF (beta == z_zero) THEN
378 :
379 0 : IF (alpha == z_zero) THEN
380 0 : a(:, :) = z_zero
381 0 : ELSE IF (alpha == z_one) THEN
382 0 : CALL timestop(handle)
383 0 : RETURN
384 : ELSE
385 0 : a(:, :) = alpha*a(:, :)
386 : END IF
387 :
388 : ELSE
389 96016 : IF (matrix_a%matrix_struct%context /= matrix_b%matrix_struct%context) &
390 0 : CPABORT("matrices must be in the same blacs context")
391 :
392 96016 : IF (cp_fm_struct_equivalent(matrix_a%matrix_struct, &
393 : matrix_b%matrix_struct)) THEN
394 :
395 96016 : b => matrix_b%local_data
396 :
397 96016 : IF (alpha == z_zero) THEN
398 7426 : IF (beta == z_one) THEN
399 : !a(:, :) = b(:, :)
400 62302 : DO icol_local = 1, ncol_local
401 960076 : DO irow_local = 1, nrow_local
402 952686 : a(irow_local, icol_local) = b(irow_local, icol_local)
403 : END DO
404 : END DO
405 : ELSE
406 : !a(:, :) = beta*b(:, :)
407 684 : DO icol_local = 1, ncol_local
408 6516 : DO irow_local = 1, nrow_local
409 6480 : a(irow_local, icol_local) = beta*b(irow_local, icol_local)
410 : END DO
411 : END DO
412 : END IF
413 88590 : ELSE IF (alpha == z_one) THEN
414 23844 : IF (beta == z_one) THEN
415 : !a(:, :) = a(:, :)+b(:, :)
416 178296 : DO icol_local = 1, ncol_local
417 2273536 : DO irow_local = 1, nrow_local
418 2264592 : a(irow_local, icol_local) = a(irow_local, icol_local) + b(irow_local, icol_local)
419 : END DO
420 : END DO
421 : ELSE
422 : !a(:, :) = a(:, :)+beta*b(:, :)
423 402432 : DO icol_local = 1, ncol_local
424 11836780 : DO irow_local = 1, nrow_local
425 11821880 : a(irow_local, icol_local) = a(irow_local, icol_local) + beta*b(irow_local, icol_local)
426 : END DO
427 : END DO
428 : END IF
429 : ELSE
430 : !a(:, :) = alpha*a(:, :)+beta*b(:, :)
431 562290 : DO icol_local = 1, ncol_local
432 5184146 : DO irow_local = 1, nrow_local
433 5119400 : a(irow_local, icol_local) = alpha*a(irow_local, icol_local) + beta*b(irow_local, icol_local)
434 : END DO
435 : END DO
436 : END IF
437 : ELSE
438 : CALL cp_abort(__LOCATION__, &
439 : "cp_cfm_scale_and_add_fm is not yet implemented for cases "// &
440 0 : "where two input matrix structures are not equivalent")
441 : END IF
442 : END IF
443 96016 : CALL timestop(handle)
444 96016 : 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) /= 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) /= 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 1260922 : 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(INOUT) :: 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 1260922 : 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 1260922 : CALL timeset(routineN, handle)
583 1260922 : a => matrix_a%local_data
584 1260922 : b => matrix_b%local_data
585 1260922 : c => matrix_c%local_data
586 :
587 1260922 : i_a = 1
588 1260922 : IF (PRESENT(a_first_row)) i_a = a_first_row
589 :
590 1260922 : j_a = 1
591 1260922 : IF (PRESENT(a_first_col)) j_a = a_first_col
592 :
593 1260922 : i_b = 1
594 1260922 : IF (PRESENT(b_first_row)) i_b = b_first_row
595 :
596 1260922 : j_b = 1
597 1260922 : IF (PRESENT(b_first_col)) j_b = b_first_col
598 :
599 1260922 : i_c = 1
600 1260922 : IF (PRESENT(c_first_row)) i_c = c_first_row
601 :
602 1260922 : j_c = 1
603 1260922 : IF (PRESENT(c_first_col)) j_c = c_first_col
604 :
605 : #if defined(__parallel)
606 12609220 : desca(:) = matrix_a%matrix_struct%descriptor(:)
607 12609220 : descb(:) = matrix_b%matrix_struct%descriptor(:)
608 12609220 : 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 1260922 : 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 1260922 : CALL timestop(handle)
622 1260922 : END SUBROUTINE cp_cfm_gemm
623 :
624 : ! **************************************************************************************************
625 : !> \brief Scales columns of the full matrix by corresponding factors.
626 : !> \param matrix_a target matrix; columns outside the selected window are unchanged
627 : !> \param scaling factors for consecutive columns starting at first_col;
628 : !> the window is truncated at the last global column of the matrix
629 : !> \param first_col first global column to scale (defaults to 1)
630 : !> \param source optional source with the same shape and distribution as matrix_a;
631 : !> if omitted, matrix_a is scaled in place
632 : !> \author Joost VandeVondele
633 : ! **************************************************************************************************
634 199318 : SUBROUTINE cp_cfm_column_scale(matrix_a, scaling, first_col, source)
635 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_a
636 : COMPLEX(kind=dp), DIMENSION(:), INTENT(IN) :: scaling
637 : INTEGER, INTENT(IN), OPTIONAL :: first_col
638 : TYPE(cp_cfm_type), INTENT(IN), OPTIONAL :: source
639 :
640 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_column_scale'
641 :
642 199318 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b
643 : INTEGER :: first, handle, i_factor, icol_local, &
644 : ncol_local, nrow_local
645 199318 : INTEGER, DIMENSION(:), POINTER :: col_indices
646 : LOGICAL :: compatible_layout
647 :
648 199318 : CALL timeset(routineN, handle)
649 :
650 199318 : a => matrix_a%local_data
651 199318 : b => a
652 199318 : IF (PRESENT(source)) THEN
653 6912 : compatible_layout = cp_fm_struct_equivalent(matrix_a%matrix_struct, source%matrix_struct)
654 6912 : CPASSERT(compatible_layout)
655 6912 : b => source%local_data
656 : END IF
657 199318 : first = 1
658 199318 : IF (PRESENT(first_col)) first = first_col
659 6912 : CPASSERT(first >= 1)
660 199318 : CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, col_indices=col_indices)
661 :
662 2708815 : DO icol_local = 1, ncol_local
663 2509497 : i_factor = col_indices(icol_local) - first + 1
664 2509497 : IF (i_factor < 1 .OR. i_factor > SIZE(scaling)) CYCLE
665 192746744 : a(1:nrow_local, icol_local) = scaling(i_factor)*b(1:nrow_local, icol_local)
666 : END DO
667 :
668 199318 : CALL timestop(handle)
669 199318 : END SUBROUTINE cp_cfm_column_scale
670 :
671 : ! **************************************************************************************************
672 : !> \brief Scales a complex matrix by a real number.
673 : !> matrix_a = alpha * matrix_b
674 : !> \param alpha scale factor
675 : !> \param matrix_a complex matrix to scale
676 : ! **************************************************************************************************
677 22213 : SUBROUTINE cp_cfm_dscale(alpha, matrix_a)
678 : REAL(kind=dp), INTENT(IN) :: alpha
679 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_a
680 :
681 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_dscale'
682 :
683 22213 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
684 : INTEGER :: handle
685 :
686 22213 : CALL timeset(routineN, handle)
687 :
688 22213 : NULLIFY (a)
689 :
690 22213 : a => matrix_a%local_data
691 :
692 66639 : CALL zdscal(SIZE(a), alpha, a(1, 1), 1)
693 :
694 22213 : CALL timestop(handle)
695 22213 : END SUBROUTINE cp_cfm_dscale
696 :
697 : ! **************************************************************************************************
698 : !> \brief Scales a complex matrix by a complex number.
699 : !> matrix_a = alpha * matrix_b
700 : !> \param alpha scale factor
701 : !> \param matrix_a complex matrix to scale
702 : !> \note
703 : !> use cp_fm_set_all to zero (avoids problems with nan)
704 : ! **************************************************************************************************
705 39645 : SUBROUTINE cp_cfm_zscale(alpha, matrix_a)
706 : COMPLEX(kind=dp), INTENT(IN) :: alpha
707 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_a
708 :
709 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_zscale'
710 :
711 39645 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
712 : INTEGER :: handle
713 :
714 39645 : CALL timeset(routineN, handle)
715 :
716 39645 : NULLIFY (a)
717 :
718 39645 : a => matrix_a%local_data
719 :
720 6460567 : a(:, :) = alpha*a(:, :)
721 :
722 39645 : CALL timestop(handle)
723 39645 : END SUBROUTINE cp_cfm_zscale
724 :
725 : ! **************************************************************************************************
726 : !> \brief Solve the system of linear equations A*b=A_general using LU decomposition.
727 : !> Pay attention that both matrices are overwritten on exit and that
728 : !> the result is stored into the matrix 'general_a'.
729 : !> \param matrix_a matrix A (overwritten on exit)
730 : !> \param general_a (input) matrix A_general, (output) matrix B
731 : !> \param determinant (optional) determinant
732 : !> \author Florian Schiffmann
733 : ! **************************************************************************************************
734 7144 : SUBROUTINE cp_cfm_solve(matrix_a, general_a, determinant)
735 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, general_a
736 : COMPLEX(kind=dp), OPTIONAL :: determinant
737 :
738 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_solve'
739 :
740 7144 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, a_general
741 : INTEGER :: counter, handle, info, irow, nrow_global
742 7144 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
743 :
744 : #if defined(__parallel)
745 : INTEGER :: icol, ncol_local, nrow_local
746 : INTEGER, DIMENSION(9) :: desca, descb
747 7144 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
748 : #else
749 : INTEGER :: lda, ldb
750 : #endif
751 :
752 7144 : CALL timeset(routineN, handle)
753 :
754 7144 : a => matrix_a%local_data
755 7144 : a_general => general_a%local_data
756 7144 : nrow_global = matrix_a%matrix_struct%nrow_global
757 21432 : ALLOCATE (ipivot(nrow_global))
758 :
759 : #if defined(__parallel)
760 71440 : desca(:) = matrix_a%matrix_struct%descriptor(:)
761 71440 : descb(:) = general_a%matrix_struct%descriptor(:)
762 7144 : CALL pzgetrf(nrow_global, nrow_global, a(1, 1), 1, 1, desca, ipivot, info)
763 7144 : IF (PRESENT(determinant)) THEN
764 : CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, &
765 6416 : row_indices=row_indices, col_indices=col_indices)
766 :
767 6416 : counter = 0
768 19296 : DO irow = 1, nrow_local
769 19296 : IF (ipivot(irow) /= row_indices(irow)) counter = counter + 1
770 : END DO
771 :
772 6416 : IF (MOD(counter, 2) == 0) THEN
773 6406 : determinant = z_one
774 : ELSE
775 10 : determinant = -z_one
776 : END IF
777 :
778 : ! compute product of diagonal elements
779 : irow = 1
780 : icol = 1
781 28932 : DO WHILE (irow <= nrow_local .AND. icol <= ncol_local)
782 28932 : IF (row_indices(irow) < col_indices(icol)) THEN
783 0 : irow = irow + 1
784 22516 : ELSE IF (row_indices(irow) > col_indices(icol)) THEN
785 9636 : icol = icol + 1
786 : ELSE ! diagonal element
787 12880 : determinant = determinant*a(irow, icol)
788 12880 : irow = irow + 1
789 12880 : icol = icol + 1
790 : END IF
791 : END DO
792 6416 : CALL matrix_a%matrix_struct%para_env%prod(determinant)
793 : END IF
794 :
795 : CALL pzgetrs("N", nrow_global, nrow_global, a(1, 1), 1, 1, desca, &
796 7144 : ipivot, a_general(1, 1), 1, 1, descb, info)
797 : #else
798 : lda = SIZE(a, 1)
799 : ldb = SIZE(a_general, 1)
800 : CALL zgetrf(nrow_global, nrow_global, a(1, 1), lda, ipivot, info)
801 : IF (PRESENT(determinant)) THEN
802 : counter = 0
803 : determinant = z_one
804 : DO irow = 1, nrow_global
805 : IF (ipivot(irow) /= irow) counter = counter + 1
806 : determinant = determinant*a(irow, irow)
807 : END DO
808 : IF (MOD(counter, 2) == 1) determinant = -1.0_dp*determinant
809 : END IF
810 : CALL zgetrs("N", nrow_global, nrow_global, a(1, 1), lda, ipivot, a_general(1, 1), ldb, info)
811 : #endif
812 :
813 : ! info is allowed to be zero
814 : ! this does just signal a zero diagonal element
815 7144 : DEALLOCATE (ipivot)
816 7144 : CALL timestop(handle)
817 :
818 7144 : END SUBROUTINE cp_cfm_solve
819 :
820 : ! **************************************************************************************************
821 : !> \brief Inverts a matrix using LU decomposition. The input matrix will be overwritten.
822 : !> \param matrix input a general square non-singular matrix, outputs its inverse
823 : !> \param info_out optional, if present outputs the info from (p)zgetri
824 : !> \author Lianheng Tong
825 : ! **************************************************************************************************
826 120050 : SUBROUTINE cp_cfm_lu_invert(matrix, info_out)
827 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
828 : INTEGER, INTENT(out), OPTIONAL :: info_out
829 :
830 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_lu_invert'
831 :
832 120050 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
833 : COMPLEX(kind=dp), DIMENSION(1) :: work1
834 120050 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: mat
835 : INTEGER :: handle, info, lwork, nrows_global
836 120050 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
837 :
838 : #if defined(__parallel)
839 : INTEGER :: liwork
840 120050 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
841 : INTEGER, DIMENSION(1) :: iwork1
842 : INTEGER, DIMENSION(9) :: desca
843 : #else
844 : INTEGER :: lda
845 : #endif
846 :
847 120050 : CALL timeset(routineN, handle)
848 :
849 120050 : mat => matrix%local_data
850 120050 : nrows_global = matrix%matrix_struct%nrow_global
851 120050 : CPASSERT(nrows_global == matrix%matrix_struct%ncol_global)
852 360150 : ALLOCATE (ipivot(nrows_global))
853 :
854 : ! do LU decomposition
855 : #if defined(__parallel)
856 1200500 : desca = matrix%matrix_struct%descriptor
857 : CALL pzgetrf(nrows_global, nrows_global, &
858 120050 : mat(1, 1), 1, 1, desca, ipivot, info)
859 : #else
860 : lda = SIZE(mat, 1)
861 : CALL zgetrf(nrows_global, nrows_global, &
862 : mat(1, 1), lda, ipivot, info)
863 : #endif
864 120050 : IF (info /= 0) THEN
865 0 : CALL cp_abort(__LOCATION__, "LU decomposition has failed")
866 : END IF
867 :
868 : ! do inversion
869 : #if defined(__parallel)
870 : CALL pzgetri(nrows_global, mat(1, 1), 1, 1, desca, &
871 120050 : ipivot, work1, -1, iwork1, -1, info)
872 120050 : lwork = INT(work1(1))
873 120050 : liwork = INT(iwork1(1))
874 360150 : ALLOCATE (work(lwork))
875 360150 : ALLOCATE (iwork(liwork))
876 : CALL pzgetri(nrows_global, mat(1, 1), 1, 1, desca, &
877 120050 : ipivot, work, lwork, iwork, liwork, info)
878 120050 : DEALLOCATE (iwork)
879 : #else
880 : CALL zgetri(nrows_global, mat(1, 1), lda, ipivot, work1, -1, info)
881 : lwork = INT(work1(1))
882 : ALLOCATE (work(lwork))
883 : CALL zgetri(nrows_global, mat(1, 1), lda, ipivot, work, lwork, info)
884 : #endif
885 120050 : DEALLOCATE (work)
886 120050 : DEALLOCATE (ipivot)
887 :
888 120050 : IF (PRESENT(info_out)) THEN
889 0 : info_out = info
890 : ELSE
891 120050 : IF (info /= 0) &
892 0 : CALL cp_abort(__LOCATION__, "LU inversion has failed")
893 : END IF
894 :
895 120050 : CALL timestop(handle)
896 :
897 120050 : END SUBROUTINE cp_cfm_lu_invert
898 :
899 : ! **************************************************************************************************
900 : !> \brief By default, returns the trace of matrix_a^T matrix_b, i.e
901 : !> sum_{i,j}(matrix_a(i,j)*matrix_b(i,j)) .
902 : !> \param matrix_a a complex matrix
903 : !> \param matrix_b another complex matrix
904 : !> \param trace value of the trace operator
905 : !> \param conjugate_a conjugate the first matrix for a Hermitian Frobenius inner product (default: false)
906 : !> \param local return only the local contribution, leaving the reduction to the caller (default: false)
907 : !> \par History
908 : !> * 09.2017 created [Sergey Chulkov]
909 : !> \author Sergey Chulkov
910 : !> \note
911 : !> Based on the subroutine cp_fm_trace(). Note the transposition of matrix_a!
912 : ! **************************************************************************************************
913 167358 : SUBROUTINE cp_cfm_trace(matrix_a, matrix_b, trace, conjugate_a, local)
914 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, matrix_b
915 : COMPLEX(kind=dp), INTENT(out) :: trace
916 : LOGICAL, INTENT(IN), OPTIONAL :: conjugate_a, local
917 :
918 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_trace'
919 :
920 : INTEGER :: handle, mypcol, myprow, ncol_local, &
921 : npcol, nprow, nrow_local
922 : LOGICAL :: reduce
923 : TYPE(cp_blacs_env_type), POINTER :: context
924 :
925 167358 : CALL timeset(routineN, handle)
926 :
927 167358 : context => matrix_a%matrix_struct%context
928 167358 : myprow = context%mepos(1)
929 167358 : mypcol = context%mepos(2)
930 167358 : nprow = context%num_pe(1)
931 167358 : npcol = context%num_pe(2)
932 :
933 167358 : nrow_local = MIN(matrix_a%matrix_struct%nrow_locals(myprow), matrix_b%matrix_struct%nrow_locals(myprow))
934 167358 : ncol_local = MIN(matrix_a%matrix_struct%ncol_locals(mypcol), matrix_b%matrix_struct%ncol_locals(mypcol))
935 :
936 : ! compute an accurate dot-product
937 : trace = accurate_dot_product(matrix_a%local_data(1:nrow_local, 1:ncol_local), &
938 167358 : matrix_b%local_data(1:nrow_local, 1:ncol_local), conjugate_a=conjugate_a)
939 :
940 167358 : reduce = .TRUE.
941 167358 : IF (PRESENT(local)) reduce = .NOT. local
942 167358 : IF (reduce) CALL matrix_a%matrix_struct%para_env%sum(trace)
943 :
944 167358 : CALL timestop(handle)
945 :
946 167358 : END SUBROUTINE cp_cfm_trace
947 :
948 : ! **************************************************************************************************
949 : !> \brief Multiplies in place by a triangular matrix:
950 : !> matrix_b = alpha op(triangular_matrix) matrix_b
951 : !> or (if side='R')
952 : !> matrix_b = alpha matrix_b op(triangular_matrix)
953 : !> op(triangular_matrix) is:
954 : !> triangular_matrix (if transa="N" and invert_tr=.false.)
955 : !> triangular_matrix^T (if transa="T" and invert_tr=.false.)
956 : !> triangular_matrix^H (if transa="C" and invert_tr=.false.)
957 : !> triangular_matrix^(-1) (if transa="N" and invert_tr=.true.)
958 : !> triangular_matrix^(-T) (if transa="T" and invert_tr=.true.)
959 : !> triangular_matrix^(-H) (if transa="C" and invert_tr=.true.)
960 : !> \param triangular_matrix the triangular matrix that multiplies the other
961 : !> \param matrix_b the matrix that gets multiplied and stores the result
962 : !> \param side on which side of matrix_b stays op(triangular_matrix)
963 : !> (defaults to 'L')
964 : !> \param transa_tr ...
965 : !> \param invert_tr if the triangular matrix should be inverted
966 : !> (defaults to false)
967 : !> \param uplo_tr if triangular_matrix is stored in the upper ('U') or
968 : !> lower ('L') triangle (defaults to 'U')
969 : !> \param unit_diag_tr if the diagonal elements of triangular_matrix should
970 : !> be assumed to be 1 (defaults to false)
971 : !> \param n_rows the number of rows of the result (defaults to
972 : !> size(matrix_b,1))
973 : !> \param n_cols the number of columns of the result (defaults to
974 : !> size(matrix_b,2))
975 : !> \param alpha ...
976 : !> \par History
977 : !> 08.2002 created [fawzi]
978 : !> \author Fawzi Mohamed
979 : !> \note
980 : !> needs an mpi env
981 : ! **************************************************************************************************
982 240866 : SUBROUTINE cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, &
983 : transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, &
984 : alpha)
985 : TYPE(cp_cfm_type), INTENT(IN) :: triangular_matrix, matrix_b
986 : CHARACTER, INTENT(in), OPTIONAL :: side, transa_tr
987 : LOGICAL, INTENT(in), OPTIONAL :: invert_tr
988 : CHARACTER, INTENT(in), OPTIONAL :: uplo_tr
989 : LOGICAL, INTENT(in), OPTIONAL :: unit_diag_tr
990 : INTEGER, INTENT(in), OPTIONAL :: n_rows, n_cols
991 : COMPLEX(kind=dp), INTENT(in), OPTIONAL :: alpha
992 :
993 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_triangular_multiply'
994 :
995 : CHARACTER :: side_char, transa, unit_diag, uplo
996 : COMPLEX(kind=dp) :: al
997 : INTEGER :: handle, m, n
998 : LOGICAL :: invert
999 :
1000 120433 : CALL timeset(routineN, handle)
1001 120433 : side_char = 'L'
1002 120433 : unit_diag = 'N'
1003 120433 : uplo = 'U'
1004 120433 : transa = 'N'
1005 120433 : invert = .FALSE.
1006 120433 : al = z_one
1007 120433 : CALL cp_cfm_get_info(matrix_b, nrow_global=m, ncol_global=n)
1008 120433 : IF (PRESENT(side)) side_char = side
1009 120433 : IF (PRESENT(invert_tr)) invert = invert_tr
1010 120433 : IF (PRESENT(uplo_tr)) uplo = uplo_tr
1011 120433 : IF (PRESENT(unit_diag_tr)) THEN
1012 0 : IF (unit_diag_tr) THEN
1013 0 : unit_diag = 'U'
1014 : ELSE
1015 : unit_diag = 'N'
1016 : END IF
1017 : END IF
1018 120433 : IF (PRESENT(transa_tr)) transa = transa_tr
1019 120433 : IF (PRESENT(alpha)) al = alpha
1020 120433 : IF (PRESENT(n_rows)) m = n_rows
1021 120433 : IF (PRESENT(n_cols)) n = n_cols
1022 :
1023 120433 : IF (invert) THEN
1024 :
1025 : #if defined(__parallel)
1026 : CALL pztrsm(side_char, uplo, transa, unit_diag, m, n, al, &
1027 : triangular_matrix%local_data(1, 1), 1, 1, &
1028 : triangular_matrix%matrix_struct%descriptor, &
1029 : matrix_b%local_data(1, 1), 1, 1, &
1030 13982 : matrix_b%matrix_struct%descriptor(1))
1031 : #else
1032 : CALL ztrsm(side_char, uplo, transa, unit_diag, m, n, al, &
1033 : triangular_matrix%local_data(1, 1), &
1034 : SIZE(triangular_matrix%local_data, 1), &
1035 : matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1))
1036 : #endif
1037 :
1038 : ELSE
1039 :
1040 : #if defined(__parallel)
1041 : CALL pztrmm(side_char, uplo, transa, unit_diag, m, n, al, &
1042 : triangular_matrix%local_data(1, 1), 1, 1, &
1043 : triangular_matrix%matrix_struct%descriptor, &
1044 : matrix_b%local_data(1, 1), 1, 1, &
1045 106451 : matrix_b%matrix_struct%descriptor(1))
1046 : #else
1047 : CALL ztrmm(side_char, uplo, transa, unit_diag, m, n, al, &
1048 : triangular_matrix%local_data(1, 1), &
1049 : SIZE(triangular_matrix%local_data, 1), &
1050 : matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1))
1051 : #endif
1052 :
1053 : END IF
1054 :
1055 120433 : CALL timestop(handle)
1056 :
1057 120433 : END SUBROUTINE cp_cfm_triangular_multiply
1058 :
1059 : ! **************************************************************************************************
1060 : !> \brief Inverts a triangular matrix.
1061 : !> \param matrix_a ...
1062 : !> \param uplo ...
1063 : !> \param info_out ...
1064 : !> \author MI
1065 : ! **************************************************************************************************
1066 35229 : SUBROUTINE cp_cfm_triangular_invert(matrix_a, uplo, info_out)
1067 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
1068 : CHARACTER, INTENT(in), OPTIONAL :: uplo
1069 : INTEGER, INTENT(out), OPTIONAL :: info_out
1070 :
1071 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_triangular_invert'
1072 :
1073 : CHARACTER :: unit_diag, my_uplo
1074 : INTEGER :: handle, info, ncol_global
1075 : COMPLEX(kind=dp), DIMENSION(:, :), &
1076 35229 : POINTER :: a
1077 : #if defined(__parallel)
1078 : INTEGER, DIMENSION(9) :: desca
1079 : #endif
1080 :
1081 35229 : CALL timeset(routineN, handle)
1082 :
1083 35229 : unit_diag = 'N'
1084 35229 : my_uplo = 'U'
1085 35229 : IF (PRESENT(uplo)) my_uplo = uplo
1086 :
1087 35229 : ncol_global = matrix_a%matrix_struct%ncol_global
1088 :
1089 35229 : a => matrix_a%local_data
1090 :
1091 : #if defined(__parallel)
1092 352290 : desca(:) = matrix_a%matrix_struct%descriptor(:)
1093 35229 : CALL pztrtri(my_uplo, unit_diag, ncol_global, a(1, 1), 1, 1, desca, info)
1094 : #else
1095 : CALL ztrtri(my_uplo, unit_diag, ncol_global, a(1, 1), ncol_global, info)
1096 : #endif
1097 :
1098 35229 : IF (PRESENT(info_out)) THEN
1099 0 : info_out = info
1100 : ELSE
1101 35229 : IF (info /= 0) &
1102 : CALL cp_abort(__LOCATION__, &
1103 0 : "triangular invert failed: matrix is not positive definite or ill-conditioned")
1104 : END IF
1105 :
1106 35229 : CALL timestop(handle)
1107 35229 : END SUBROUTINE cp_cfm_triangular_invert
1108 :
1109 : ! **************************************************************************************************
1110 : !> \brief Transposes a BLACS distributed complex matrix.
1111 : !> \param matrix input matrix
1112 : !> \param trans 'T' for transpose, 'C' for Hermitian conjugate
1113 : !> \param matrixt output matrix
1114 : !> \author Lianheng Tong
1115 : ! **************************************************************************************************
1116 33393 : SUBROUTINE cp_cfm_transpose(matrix, trans, matrixt)
1117 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1118 : CHARACTER, INTENT(in) :: trans
1119 : TYPE(cp_cfm_type), INTENT(IN) :: matrixt
1120 :
1121 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_transpose'
1122 :
1123 33393 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa, cc
1124 : INTEGER :: handle, ncol_global, nrow_global
1125 : #if defined(__parallel)
1126 : INTEGER, DIMENSION(9) :: desca, descc
1127 : #elif !defined(__MKL)
1128 : INTEGER :: ii, jj
1129 : #endif
1130 :
1131 33393 : CALL timeset(routineN, handle)
1132 :
1133 33393 : nrow_global = matrix%matrix_struct%nrow_global
1134 33393 : ncol_global = matrix%matrix_struct%ncol_global
1135 :
1136 33393 : CPASSERT(matrixt%matrix_struct%nrow_global == ncol_global)
1137 33393 : CPASSERT(matrixt%matrix_struct%ncol_global == nrow_global)
1138 :
1139 33393 : aa => matrix%local_data
1140 33393 : cc => matrixt%local_data
1141 :
1142 : #if defined(__parallel)
1143 333930 : desca = matrix%matrix_struct%descriptor
1144 333930 : descc = matrixt%matrix_struct%descriptor
1145 13338 : SELECT CASE (trans)
1146 : CASE ('T')
1147 : CALL pztranu(nrow_global, ncol_global, &
1148 : z_one, aa(1, 1), 1, 1, desca, &
1149 13338 : z_zero, cc(1, 1), 1, 1, descc)
1150 : CASE ('C')
1151 : CALL pztranc(nrow_global, ncol_global, &
1152 : z_one, aa(1, 1), 1, 1, desca, &
1153 20055 : z_zero, cc(1, 1), 1, 1, descc)
1154 : CASE DEFAULT
1155 33393 : CPABORT("trans only accepts 'T' or 'C'")
1156 : END SELECT
1157 : #elif defined(__MKL)
1158 : CALL mkl_zomatcopy('C', trans, nrow_global, ncol_global, 1.0_dp, aa(1, 1), nrow_global, cc(1, 1), ncol_global)
1159 : #else
1160 : SELECT CASE (trans)
1161 : CASE ('T')
1162 : DO jj = 1, ncol_global
1163 : DO ii = 1, nrow_global
1164 : cc(ii, jj) = aa(jj, ii)
1165 : END DO
1166 : END DO
1167 : CASE ('C')
1168 : DO jj = 1, ncol_global
1169 : DO ii = 1, nrow_global
1170 : cc(ii, jj) = CONJG(aa(jj, ii))
1171 : END DO
1172 : END DO
1173 : CASE DEFAULT
1174 : CPABORT("trans only accepts 'T' or 'C'")
1175 : END SELECT
1176 : #endif
1177 :
1178 33393 : CALL timestop(handle)
1179 33393 : END SUBROUTINE cp_cfm_transpose
1180 :
1181 : ! **************************************************************************************************
1182 : !> \brief Norm of matrix using (p)zlange.
1183 : !> \param matrix input a general matrix
1184 : !> \param mode 'M' max abs element value,
1185 : !> '1' or 'O' one norm, i.e. maximum column sum,
1186 : !> 'I' infinity norm, i.e. maximum row sum,
1187 : !> 'F' or 'E' Frobenius norm, i.e. sqrt of sum of all squares of elements
1188 : !> \return the norm according to mode
1189 : !> \author Lianheng Tong
1190 : ! **************************************************************************************************
1191 237266 : FUNCTION cp_cfm_norm(matrix, mode) RESULT(res)
1192 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1193 : CHARACTER, INTENT(IN) :: mode
1194 : REAL(kind=dp) :: res
1195 :
1196 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_norm'
1197 :
1198 237266 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa
1199 : INTEGER :: handle, lwork, ncols, ncols_local, &
1200 : nrows, nrows_local
1201 237266 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
1202 :
1203 : #if defined(__parallel)
1204 : INTEGER, DIMENSION(9) :: desca
1205 : #else
1206 : INTEGER :: lda
1207 : #endif
1208 :
1209 237266 : CALL timeset(routineN, handle)
1210 :
1211 : CALL cp_cfm_get_info(matrix=matrix, &
1212 : nrow_global=nrows, &
1213 : ncol_global=ncols, &
1214 : nrow_local=nrows_local, &
1215 237266 : ncol_local=ncols_local)
1216 237266 : aa => matrix%local_data
1217 :
1218 : SELECT CASE (mode)
1219 : CASE ('M', 'm')
1220 0 : lwork = 1
1221 : CASE ('1', 'O', 'o')
1222 : #if defined(__parallel)
1223 0 : lwork = ncols_local
1224 : #else
1225 : lwork = 1
1226 : #endif
1227 : CASE ('I', 'i')
1228 : #if defined(__parallel)
1229 8 : lwork = nrows_local
1230 : #else
1231 : lwork = nrows
1232 : #endif
1233 : CASE ('F', 'f', 'E', 'e')
1234 0 : lwork = 1
1235 : CASE DEFAULT
1236 237266 : CPABORT("mode input is not valid")
1237 : END SELECT
1238 :
1239 711798 : ALLOCATE (work(lwork))
1240 :
1241 : #if defined(__parallel)
1242 2372660 : desca = matrix%matrix_struct%descriptor
1243 237266 : res = pzlange(mode, nrows, ncols, aa(1, 1), 1, 1, desca, work)
1244 : #else
1245 : lda = SIZE(aa, 1)
1246 : res = zlange(mode, nrows, ncols, aa(1, 1), lda, work)
1247 : #endif
1248 :
1249 237266 : DEALLOCATE (work)
1250 237266 : CALL timestop(handle)
1251 237266 : END FUNCTION cp_cfm_norm
1252 :
1253 : ! **************************************************************************************************
1254 : !> \brief Applies a planar rotation defined by cs and sn to the i'th and j'th rows.
1255 : !> \param matrix ...
1256 : !> \param irow ...
1257 : !> \param jrow ...
1258 : !> \param cs cosine of the rotation angle
1259 : !> \param sn sinus of the rotation angle
1260 : !> \author Ole Schuett
1261 : ! **************************************************************************************************
1262 375120 : SUBROUTINE cp_cfm_rot_rows(matrix, irow, jrow, cs, sn)
1263 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1264 : INTEGER, INTENT(IN) :: irow, jrow
1265 : REAL(dp), INTENT(IN) :: cs, sn
1266 :
1267 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_rot_rows'
1268 : INTEGER :: handle, ncol
1269 : COMPLEX(KIND=dp) :: sn_cmplx
1270 :
1271 : #if defined(__parallel)
1272 : INTEGER :: info, lwork
1273 : INTEGER, DIMENSION(9) :: desc
1274 375120 : REAL(dp), DIMENSION(:), ALLOCATABLE :: work
1275 : #endif
1276 375120 : CALL timeset(routineN, handle)
1277 375120 : CALL cp_cfm_get_info(matrix, ncol_global=ncol)
1278 375120 : sn_cmplx = CMPLX(sn, 0.0_dp, dp)
1279 : #if defined(__parallel)
1280 375120 : IF (1 /= matrix%matrix_struct%context%n_pid) THEN
1281 375120 : lwork = 2*ncol + 1
1282 1125360 : ALLOCATE (work(lwork))
1283 3751200 : desc(:) = matrix%matrix_struct%descriptor(:)
1284 375120 : info = 0
1285 : CALL pzrot(ncol, &
1286 : matrix%local_data(1, 1), irow, 1, desc, ncol, &
1287 : matrix%local_data(1, 1), jrow, 1, desc, ncol, &
1288 375120 : cs, sn_cmplx, work, lwork, info)
1289 375120 : CPASSERT(info == 0)
1290 375120 : DEALLOCATE (work)
1291 : ELSE
1292 : #endif
1293 0 : CALL zrot(ncol, matrix%local_data(irow, 1), ncol, matrix%local_data(jrow, 1), ncol, cs, sn_cmplx)
1294 : #if defined(__parallel)
1295 : END IF
1296 : #endif
1297 375120 : CALL timestop(handle)
1298 375120 : END SUBROUTINE cp_cfm_rot_rows
1299 :
1300 : ! **************************************************************************************************
1301 : !> \brief Applies a planar rotation defined by cs and sn to the i'th and j'th columnns.
1302 : !> \param matrix ...
1303 : !> \param icol ...
1304 : !> \param jcol ...
1305 : !> \param cs cosine of the rotation angle
1306 : !> \param sn sinus of the rotation angle
1307 : !> \author Ole Schuett
1308 : ! **************************************************************************************************
1309 422760 : SUBROUTINE cp_cfm_rot_cols(matrix, icol, jcol, cs, sn)
1310 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1311 : INTEGER, INTENT(IN) :: icol, jcol
1312 : REAL(dp), INTENT(IN) :: cs, sn
1313 :
1314 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_rot_cols'
1315 : INTEGER :: handle, nrow
1316 : COMPLEX(KIND=dp) :: sn_cmplx
1317 :
1318 : #if defined(__parallel)
1319 : INTEGER :: info, lwork
1320 : INTEGER, DIMENSION(9) :: desc
1321 422760 : REAL(dp), DIMENSION(:), ALLOCATABLE :: work
1322 : #endif
1323 422760 : CALL timeset(routineN, handle)
1324 422760 : CALL cp_cfm_get_info(matrix, nrow_global=nrow)
1325 422760 : sn_cmplx = CMPLX(sn, 0.0_dp, dp)
1326 : #if defined(__parallel)
1327 422760 : IF (1 /= matrix%matrix_struct%context%n_pid) THEN
1328 422760 : lwork = 2*nrow + 1
1329 1268280 : ALLOCATE (work(lwork))
1330 4227600 : desc(:) = matrix%matrix_struct%descriptor(:)
1331 422760 : info = 0
1332 : CALL pzrot(nrow, &
1333 : matrix%local_data(1, 1), 1, icol, desc, 1, &
1334 : matrix%local_data(1, 1), 1, jcol, desc, 1, &
1335 422760 : cs, sn_cmplx, work, lwork, info)
1336 422760 : CPASSERT(info == 0)
1337 422760 : DEALLOCATE (work)
1338 : ELSE
1339 : #endif
1340 0 : CALL zrot(nrow, matrix%local_data(1, icol), 1, matrix%local_data(1, jcol), 1, cs, sn_cmplx)
1341 : #if defined(__parallel)
1342 : END IF
1343 : #endif
1344 422760 : CALL timestop(handle)
1345 422760 : END SUBROUTINE cp_cfm_rot_cols
1346 :
1347 : ! **************************************************************************************************
1348 : !> \brief ...
1349 : !> \param matrix ...
1350 : !> \param workspace ...
1351 : !> \param uplo triangular format; defaults to 'U'
1352 : !> \par History
1353 : !> 12.2024 Added optional workspace as input [Rocco Meli]
1354 : !> \author Jan Wilhelm
1355 : ! **************************************************************************************************
1356 15122 : SUBROUTINE cp_cfm_uplo_to_full(matrix, workspace, uplo)
1357 :
1358 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1359 : TYPE(cp_cfm_type), INTENT(IN), OPTIONAL :: workspace
1360 : CHARACTER, INTENT(IN), OPTIONAL :: uplo
1361 :
1362 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_cfm_uplo_to_full'
1363 :
1364 : CHARACTER :: myuplo
1365 : INTEGER :: handle, i_global, iiB, j_global, jjB, &
1366 : ncol_local, nrow_local
1367 7561 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1368 : TYPE(cp_cfm_type) :: work
1369 :
1370 7561 : CALL timeset(routineN, handle)
1371 :
1372 7561 : IF (.NOT. PRESENT(workspace)) THEN
1373 6910 : CALL cp_cfm_create(work, matrix%matrix_struct)
1374 : ELSE
1375 651 : work = workspace
1376 : END IF
1377 :
1378 7561 : myuplo = 'U'
1379 7561 : IF (PRESENT(uplo)) myuplo = uplo
1380 :
1381 : ! get info of fm_mat_Q
1382 : CALL cp_cfm_get_info(matrix=matrix, &
1383 : nrow_local=nrow_local, &
1384 : ncol_local=ncol_local, &
1385 : row_indices=row_indices, &
1386 7561 : col_indices=col_indices)
1387 :
1388 278263 : DO jjB = 1, ncol_local
1389 270702 : j_global = col_indices(jjB)
1390 9086380 : DO iiB = 1, nrow_local
1391 8808117 : i_global = row_indices(iiB)
1392 9078819 : IF (MERGE(j_global < i_global, j_global > i_global, (myuplo == "U") .OR. (myuplo == "u"))) THEN
1393 4331539 : matrix%local_data(iiB, jjB) = z_zero
1394 4476578 : ELSE IF (j_global == i_global) THEN
1395 145039 : matrix%local_data(iiB, jjB) = matrix%local_data(iiB, jjB)/(2.0_dp, 0.0_dp)
1396 : END IF
1397 : END DO
1398 : END DO
1399 :
1400 7561 : CALL cp_cfm_transpose(matrix, 'C', work)
1401 :
1402 7561 : CALL cp_cfm_scale_and_add(z_one, matrix, z_one, work)
1403 :
1404 7561 : IF (.NOT. PRESENT(workspace)) THEN
1405 6910 : CALL cp_cfm_release(work)
1406 : END IF
1407 :
1408 7561 : CALL timestop(handle)
1409 :
1410 7561 : END SUBROUTINE cp_cfm_uplo_to_full
1411 :
1412 : ! **************************************************************************************************
1413 : !> \brief find the norm of each column norm_{j}= sqrt( \sum_{i} A_{ij}*conjg(A_{ij}) )
1414 : !> Complex-valued mirror of cp_fm_vectorsnorm.
1415 : !> \param matrix ...
1416 : !> \param norm_array ...
1417 : !> \param squared return squared norms instead (default: false)
1418 : ! **************************************************************************************************
1419 190410 : SUBROUTINE cp_cfm_vectorsnorm(matrix, norm_array, squared)
1420 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1421 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: norm_array
1422 : LOGICAL, INTENT(IN), OPTIONAL :: squared
1423 :
1424 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_cfm_vectorsnorm'
1425 :
1426 : INTEGER :: handle, i, j, ncol_local, nrow_local
1427 :
1428 : !$ INTEGER :: nworkers
1429 190410 : INTEGER, DIMENSION(:), POINTER :: col_indices
1430 : LOGICAL :: take_sqrt
1431 : REAL(KIND=dp) :: norm_sq
1432 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1433 190410 : POINTER :: local_data
1434 :
1435 190410 : CALL timeset(routineN, handle)
1436 :
1437 190410 : local_data => matrix%local_data
1438 : CALL cp_cfm_get_info(matrix, col_indices=col_indices, nrow_local=nrow_local, &
1439 190410 : ncol_local=ncol_local)
1440 :
1441 993886 : norm_array = 0.0_dp
1442 : ! Own whole columns and keep at least 8192 elements per worker to amortize team startup.
1443 190410 : !$ nworkers = MAX(1, MIN(omp_get_max_threads(), ncol_local))
1444 190410 : !$ nworkers = MAX(1, MIN(nworkers, (nrow_local*ncol_local + 1)/8192))
1445 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) &
1446 : !$OMP SHARED(local_data, col_indices, nrow_local, ncol_local, norm_array) PRIVATE(i, j, norm_sq) &
1447 190410 : !$OMP NUM_THREADS(nworkers) IF(nworkers > 1)
1448 : DO j = 1, ncol_local
1449 : norm_sq = 0.0_dp
1450 : !$OMP SIMD REDUCTION(+:norm_sq)
1451 : DO i = 1, nrow_local
1452 : norm_sq = norm_sq + REAL(local_data(i, j), KIND=dp)**2 + AIMAG(local_data(i, j))**2
1453 : END DO
1454 : !$OMP END SIMD
1455 : norm_array(col_indices(j)) = norm_sq
1456 : END DO
1457 : !$OMP END PARALLEL DO
1458 1797362 : CALL matrix%matrix_struct%para_env%sum(norm_array)
1459 190410 : take_sqrt = .TRUE.
1460 190410 : IF (PRESENT(squared)) take_sqrt = .NOT. squared
1461 991692 : IF (take_sqrt) norm_array = SQRT(norm_array)
1462 :
1463 190410 : CALL timestop(handle)
1464 190410 : END SUBROUTINE cp_cfm_vectorsnorm
1465 :
1466 : ! **************************************************************************************************
1467 : !> \brief Adds a scalar to the diagonal of a distributed complex full matrix.
1468 : !> \param matrix input/output matrix
1469 : !> \param alpha scalar added to each diagonal element
1470 : !> \param n_active optional number of physical rows/columns
1471 : ! **************************************************************************************************
1472 23512 : SUBROUTINE cp_cfm_add_on_diag(matrix, alpha, n_active)
1473 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix
1474 : COMPLEX(KIND=dp), INTENT(IN) :: alpha
1475 : INTEGER, INTENT(IN), OPTIONAL :: n_active
1476 :
1477 : INTEGER :: i_global, n_diag, ncol_global, nrow_global
1478 : #if defined(__parallel)
1479 : INTEGER :: icol_local, ipcol, iprow, irow, irow_local, &
1480 : mypcol, myprow, npcol, nprow, nrow_local
1481 : INTEGER, DIMENSION(9) :: desca
1482 11756 : INTEGER, DIMENSION(:), POINTER :: row_indices
1483 : TYPE(cp_blacs_env_type), POINTER :: context
1484 : #endif
1485 :
1486 11756 : CALL cp_cfm_get_info(matrix, nrow_global=nrow_global, ncol_global=ncol_global)
1487 :
1488 11756 : n_diag = MIN(nrow_global, ncol_global)
1489 11756 : IF (PRESENT(n_active)) n_diag = n_active
1490 11756 : CPASSERT(n_diag >= 0 .AND. n_diag <= MIN(nrow_global, ncol_global))
1491 :
1492 : #if defined(__parallel)
1493 11756 : CALL cp_cfm_get_info(matrix, nrow_local=nrow_local, row_indices=row_indices)
1494 11756 : context => matrix%matrix_struct%context
1495 11756 : myprow = context%mepos(1)
1496 11756 : mypcol = context%mepos(2)
1497 11756 : nprow = context%num_pe(1)
1498 11756 : npcol = context%num_pe(2)
1499 117560 : desca(:) = matrix%matrix_struct%descriptor(:)
1500 :
1501 258286 : DO irow = 1, nrow_local
1502 246530 : i_global = row_indices(irow)
1503 246530 : IF (i_global > n_diag) EXIT
1504 : CALL infog2l(i_global, i_global, desca, nprow, npcol, myprow, mypcol, &
1505 246530 : irow_local, icol_local, iprow, ipcol)
1506 258286 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
1507 246530 : matrix%local_data(irow_local, icol_local) = matrix%local_data(irow_local, icol_local) + alpha
1508 : END IF
1509 : END DO
1510 : #else
1511 : DO i_global = 1, n_diag
1512 : matrix%local_data(i_global, i_global) = matrix%local_data(i_global, i_global) + alpha
1513 : END DO
1514 : #endif
1515 :
1516 11756 : END SUBROUTINE cp_cfm_add_on_diag
1517 :
1518 : ! **************************************************************************************************
1519 : !> \brief returns the diagonal of a complex full matrix: diag(i)= A_{ii}.
1520 : !> Each diagonal entry is owned by one process. The sum over the
1521 : !> process grid collects the entries.
1522 : !> \param matrix ...
1523 : !> \param diag ...
1524 : !> \param n_active optional number of leading diagonal elements to return
1525 : ! **************************************************************************************************
1526 28716 : SUBROUTINE cp_cfm_get_diag(matrix, diag, n_active)
1527 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
1528 : COMPLEX(KIND=dp), DIMENSION(:), INTENT(OUT) :: diag
1529 : INTEGER, INTENT(IN), OPTIONAL :: n_active
1530 :
1531 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_cfm_get_diag'
1532 :
1533 : INTEGER :: handle, i_global, n_diag, ncol_global, &
1534 : nrow_global
1535 : #if defined(__parallel)
1536 : INTEGER :: icol_local, ipcol, iprow, irow, irow_local, &
1537 : mypcol, myprow, npcol, nprow, nrow_local
1538 : INTEGER, DIMENSION(9) :: desca
1539 28716 : INTEGER, DIMENSION(:), POINTER :: row_indices
1540 : TYPE(cp_blacs_env_type), POINTER :: context
1541 : #endif
1542 :
1543 28716 : CALL timeset(routineN, handle)
1544 :
1545 28716 : CALL cp_cfm_get_info(matrix, nrow_global=nrow_global, ncol_global=ncol_global)
1546 :
1547 28716 : n_diag = MIN(nrow_global, ncol_global)
1548 28716 : IF (PRESENT(n_active)) n_diag = n_active
1549 28716 : CPASSERT(n_diag >= 0 .AND. n_diag <= MIN(nrow_global, ncol_global))
1550 28716 : CPASSERT(SIZE(diag) >= n_diag)
1551 :
1552 715063 : diag = z_zero
1553 : #if defined(__parallel)
1554 28716 : CALL cp_cfm_get_info(matrix, nrow_local=nrow_local, row_indices=row_indices)
1555 28716 : context => matrix%matrix_struct%context
1556 28716 : myprow = context%mepos(1)
1557 28716 : mypcol = context%mepos(2)
1558 28716 : nprow = context%num_pe(1)
1559 28716 : npcol = context%num_pe(2)
1560 287160 : desca(:) = matrix%matrix_struct%descriptor(:)
1561 :
1562 395124 : DO irow = 1, nrow_local
1563 366408 : i_global = row_indices(irow)
1564 366408 : IF (i_global > n_diag) EXIT
1565 : CALL infog2l(i_global, i_global, desca, nprow, npcol, myprow, mypcol, &
1566 366408 : irow_local, icol_local, iprow, ipcol)
1567 395124 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
1568 366408 : diag(i_global) = matrix%local_data(irow_local, icol_local)
1569 : END IF
1570 : END DO
1571 : #else
1572 : DO i_global = 1, n_diag
1573 : diag(i_global) = matrix%local_data(i_global, i_global)
1574 : END DO
1575 : #endif
1576 1401410 : CALL matrix%matrix_struct%para_env%sum(diag)
1577 :
1578 28716 : CALL timestop(handle)
1579 28716 : END SUBROUTINE cp_cfm_get_diag
1580 :
1581 : END MODULE cp_cfm_basic_linalg
|
