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