(git:f2099e5)
Loading...
Searching...
No Matches
cp_cfm_basic_linalg.F
Go to the documentation of this file.
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! **************************************************************************************************
16MODULE 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
69CONTAINS
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 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 COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a
84 INTEGER :: n, i, info, p
85 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
86 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 name="A_lu"//trim(adjustl(cp_to_string(1)))//"MATRIX")
97 CALL cp_cfm_to_cfm(matrix_a, matrix_lu)
98
99 a => matrix_lu%local_data
100 n = matrix_lu%matrix_struct%nrow_global
101 ALLOCATE (ipivot(n))
102 ipivot(:) = 0
103 p = 0
104 ALLOCATE (diag(n))
105 diag(:) = 0.0_dp
106#if defined(__parallel)
107 ! Use LU decomposition
108 desca(:) = matrix_lu%matrix_struct%descriptor(:)
109 CALL pzgetrf(n, n, a(1, 1), 1, 1, desca, ipivot, info)
110 myprow = matrix_lu%matrix_struct%context%mepos(1)
111 mypcol = matrix_lu%matrix_struct%context%mepos(2)
112 nprow = matrix_lu%matrix_struct%context%num_pe(1)
113 npcol = matrix_lu%matrix_struct%context%num_pe(2)
114 nrow_local = matrix_lu%matrix_struct%nrow_locals(myprow)
115 ncol_local = matrix_lu%matrix_struct%ncol_locals(mypcol)
116
117 DO irow_local = 1, nrow_local
118 i = matrix_lu%matrix_struct%row_indices(irow_local)
119 DO icol_local = 1, ncol_local
120 j = matrix_lu%matrix_struct%col_indices(icol_local)
121 IF (i == j) diag(i) = matrix_lu%local_data(irow_local, icol_local)
122 END DO
123 END DO
124 CALL matrix_lu%matrix_struct%para_env%sum(diag)
125 determinant = product(diag)
126 DO irow_local = 1, nrow_local
127 i = matrix_lu%matrix_struct%row_indices(irow_local)
128 IF (ipivot(irow_local) /= i) p = p + 1
129 END DO
130 CALL matrix_lu%matrix_struct%para_env%sum(p)
131 ! very important fix
132 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 DEALLOCATE (ipivot)
144 DEALLOCATE (diag)
145 CALL cp_cfm_release(matrix_lu)
146 det_a = determinant*(-2*mod(p, 2) + 1.0_dp)
147 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 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 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b, c
162 INTEGER :: handle, icol_local, irow_local, mypcol, &
163 myprow, ncol_local, nrow_local
164
165 CALL timeset(routinen, handle)
166
167 myprow = matrix_a%matrix_struct%context%mepos(1)
168 mypcol = matrix_a%matrix_struct%context%mepos(2)
169
170 a => matrix_a%local_data
171 b => matrix_b%local_data
172 c => matrix_c%local_data
173
174 nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
175 ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
176
177 DO icol_local = 1, ncol_local
178 DO irow_local = 1, nrow_local
179 c(irow_local, icol_local) = a(irow_local, icol_local)*b(irow_local, icol_local)
180 END DO
181 END DO
182
183 CALL timestop(handle)
184
185 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 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 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b, c
200 INTEGER :: handle, icol_local, irow_local, mypcol, &
201 myprow, ncol_local, nrow_local
202
203 CALL timeset(routinen, handle)
204
205 myprow = matrix_a%matrix_struct%context%mepos(1)
206 mypcol = matrix_a%matrix_struct%context%mepos(2)
207
208 a => matrix_a%local_data
209 b => matrix_b%local_data
210 c => matrix_c%local_data
211
212 nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
213 ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
214
215 DO icol_local = 1, ncol_local
216 DO irow_local = 1, nrow_local
217 c(irow_local, icol_local) = a(irow_local, icol_local)*conjg(b(irow_local, icol_local))
218 END DO
219 END DO
220
221 CALL timestop(handle)
222
223 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 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 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b
252 INTEGER :: handle, icol_local, irow_local, mypcol, &
253 myprow, ncol_local, nrow_local
254
255 CALL timeset(routinen, handle)
256
257 my_beta = z_zero
258 IF (PRESENT(beta)) my_beta = beta
259 NULLIFY (a, b)
260
261 ! to do: use dscal,dcopy,daxp
262 myprow = matrix_a%matrix_struct%context%mepos(1)
263 mypcol = matrix_a%matrix_struct%context%mepos(2)
264
265 nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
266 ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
267
268 a => matrix_a%local_data
269
270 IF (my_beta == z_zero) THEN
271
272 IF (alpha == z_zero) THEN
273 a(:, :) = z_zero
274 ELSE IF (alpha == z_one) THEN
275 CALL timestop(handle)
276 RETURN
277 ELSE
278 a(:, :) = alpha*a(:, :)
279 END IF
280
281 ELSE
282 cpassert(PRESENT(matrix_b))
283 IF (matrix_a%matrix_struct%context /= matrix_b%matrix_struct%context) &
284 cpabort("matrixes must be in the same blacs context")
285
286 IF (cp_fm_struct_equivalent(matrix_a%matrix_struct, &
287 matrix_b%matrix_struct)) THEN
288
289 b => matrix_b%local_data
290
291 IF (alpha == z_zero) THEN
292 IF (my_beta == z_one) THEN
293 !a(:, :) = b(:, :)
294 DO icol_local = 1, ncol_local
295 DO irow_local = 1, nrow_local
296 a(irow_local, icol_local) = b(irow_local, icol_local)
297 END DO
298 END DO
299 ELSE
300 !a(:, :) = my_beta*b(:, :)
301 DO icol_local = 1, ncol_local
302 DO irow_local = 1, nrow_local
303 a(irow_local, icol_local) = my_beta*b(irow_local, icol_local)
304 END DO
305 END DO
306 END IF
307 ELSE IF (alpha == z_one) THEN
308 IF (my_beta == z_one) THEN
309 !a(:, :) = a(:, :)+b(:, :)
310 DO icol_local = 1, ncol_local
311 DO irow_local = 1, nrow_local
312 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 DO icol_local = 1, ncol_local
318 DO irow_local = 1, nrow_local
319 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 DO icol_local = 1, ncol_local
326 DO irow_local = 1, nrow_local
327 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 "where input two matrix structures are not equivalent")
335 END IF
336 END IF
337 CALL timestop(handle)
338 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 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 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
360 INTEGER :: handle, icol_local, irow_local, mypcol, &
361 myprow, ncol_local, nrow_local
362 REAL(kind=dp), DIMENSION(:, :), POINTER :: b
363
364 CALL timeset(routinen, handle)
365
366 NULLIFY (a, b)
367
368 myprow = matrix_a%matrix_struct%context%mepos(1)
369 mypcol = matrix_a%matrix_struct%context%mepos(2)
370
371 nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
372 ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
373
374 a => matrix_a%local_data
375
376 IF (beta == z_zero) THEN
377
378 IF (alpha == z_zero) THEN
379 a(:, :) = z_zero
380 ELSE IF (alpha == z_one) THEN
381 CALL timestop(handle)
382 RETURN
383 ELSE
384 a(:, :) = alpha*a(:, :)
385 END IF
386
387 ELSE
388 IF (matrix_a%matrix_struct%context /= matrix_b%matrix_struct%context) &
389 cpabort("matrices must be in the same blacs context")
390
391 IF (cp_fm_struct_equivalent(matrix_a%matrix_struct, &
392 matrix_b%matrix_struct)) THEN
393
394 b => matrix_b%local_data
395
396 IF (alpha == z_zero) THEN
397 IF (beta == z_one) THEN
398 !a(:, :) = b(:, :)
399 DO icol_local = 1, ncol_local
400 DO irow_local = 1, nrow_local
401 a(irow_local, icol_local) = b(irow_local, icol_local)
402 END DO
403 END DO
404 ELSE
405 !a(:, :) = beta*b(:, :)
406 DO icol_local = 1, ncol_local
407 DO irow_local = 1, nrow_local
408 a(irow_local, icol_local) = beta*b(irow_local, icol_local)
409 END DO
410 END DO
411 END IF
412 ELSE IF (alpha == z_one) THEN
413 IF (beta == z_one) THEN
414 !a(:, :) = a(:, :)+b(:, :)
415 DO icol_local = 1, ncol_local
416 DO irow_local = 1, nrow_local
417 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 DO icol_local = 1, ncol_local
423 DO irow_local = 1, nrow_local
424 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 DO icol_local = 1, ncol_local
431 DO irow_local = 1, nrow_local
432 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 "where two input matrix structures are not equivalent")
440 END IF
441 END IF
442 CALL timestop(handle)
443 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 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 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
463 INTEGER :: counter, handle, info, irow, nrow_global
464 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
465
466#if defined(__parallel)
467 INTEGER :: icol, ncol_local, nrow_local
468 INTEGER, DIMENSION(9) :: desca
469 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
470#else
471 INTEGER :: lda
472#endif
473
474 CALL timeset(routinen, handle)
475
476 nrow_global = matrix_a%matrix_struct%nrow_global
477 a => matrix_a%local_data
478
479 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 row_indices=row_indices, col_indices=col_indices)
483
484 desca(:) = matrix_a%matrix_struct%descriptor(:)
485 CALL pzgetrf(nrow_global, nrow_global, a(1, 1), 1, 1, desca, ipivot, info)
486
487 counter = 0
488 DO irow = 1, nrow_local
489 IF (ipivot(irow) /= row_indices(irow)) counter = counter + 1
490 END DO
491
492 IF (mod(counter, 2) == 0) THEN
493 determinant = z_one
494 ELSE
495 determinant = -z_one
496 END IF
497
498 ! compute product of diagonal elements
499 irow = 1
500 icol = 1
501 DO WHILE (irow <= nrow_local .AND. icol <= ncol_local)
502 IF (row_indices(irow) < col_indices(icol)) THEN
503 irow = irow + 1
504 ELSE IF (row_indices(irow) > col_indices(icol)) THEN
505 icol = icol + 1
506 ELSE ! diagonal element
507 determinant = determinant*a(irow, icol)
508 irow = irow + 1
509 icol = icol + 1
510 END IF
511 END DO
512 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 DEALLOCATE (ipivot)
528
529 CALL timestop(handle)
530 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 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 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 CALL timeset(routinen, handle)
582 a => matrix_a%local_data
583 b => matrix_b%local_data
584 c => matrix_c%local_data
585
586 i_a = 1
587 IF (PRESENT(a_first_row)) i_a = a_first_row
588
589 j_a = 1
590 IF (PRESENT(a_first_col)) j_a = a_first_col
591
592 i_b = 1
593 IF (PRESENT(b_first_row)) i_b = b_first_row
594
595 j_b = 1
596 IF (PRESENT(b_first_col)) j_b = b_first_col
597
598 i_c = 1
599 IF (PRESENT(c_first_row)) i_c = c_first_row
600
601 j_c = 1
602 IF (PRESENT(c_first_col)) j_c = c_first_col
603
604#if defined(__parallel)
605 desca(:) = matrix_a%matrix_struct%descriptor(:)
606 descb(:) = matrix_b%matrix_struct%descriptor(:)
607 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 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 CALL timestop(handle)
621 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 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 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
638 INTEGER :: handle, icol_local, ncol_local, &
639 nrow_local
640#if defined(__parallel)
641 INTEGER, DIMENSION(:), POINTER :: col_indices
642#endif
643
644 CALL timeset(routinen, handle)
645
646 a => matrix_a%local_data
647
648#if defined(__parallel)
649 CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, col_indices=col_indices)
650 ncol_local = min(ncol_local, SIZE(scaling))
651
652 DO icol_local = 1, ncol_local
653 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 CALL timestop(handle)
665 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 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 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
680 INTEGER :: handle
681
682 CALL timeset(routinen, handle)
683
684 NULLIFY (a)
685
686 a => matrix_a%local_data
687
688 CALL zdscal(SIZE(a), alpha, a(1, 1), 1)
689
690 CALL timestop(handle)
691 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 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 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
708 INTEGER :: handle
709
710 CALL timeset(routinen, handle)
711
712 NULLIFY (a)
713
714 a => matrix_a%local_data
715
716 a(:, :) = alpha*a(:, :)
717
718 CALL timestop(handle)
719 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 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 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, a_general
737 INTEGER :: counter, handle, info, irow, nrow_global
738 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
739
740#if defined(__parallel)
741 INTEGER :: icol, ncol_local, nrow_local
742 INTEGER, DIMENSION(9) :: desca, descb
743 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
744#else
745 INTEGER :: lda, ldb
746#endif
747
748 CALL timeset(routinen, handle)
749
750 a => matrix_a%local_data
751 a_general => general_a%local_data
752 nrow_global = matrix_a%matrix_struct%nrow_global
753 ALLOCATE (ipivot(nrow_global))
754
755#if defined(__parallel)
756 desca(:) = matrix_a%matrix_struct%descriptor(:)
757 descb(:) = general_a%matrix_struct%descriptor(:)
758 CALL pzgetrf(nrow_global, nrow_global, a(1, 1), 1, 1, desca, ipivot, info)
759 IF (PRESENT(determinant)) THEN
760 CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, &
761 row_indices=row_indices, col_indices=col_indices)
762
763 counter = 0
764 DO irow = 1, nrow_local
765 IF (ipivot(irow) /= row_indices(irow)) counter = counter + 1
766 END DO
767
768 IF (mod(counter, 2) == 0) THEN
769 determinant = z_one
770 ELSE
771 determinant = -z_one
772 END IF
773
774 ! compute product of diagonal elements
775 irow = 1
776 icol = 1
777 DO WHILE (irow <= nrow_local .AND. icol <= ncol_local)
778 IF (row_indices(irow) < col_indices(icol)) THEN
779 irow = irow + 1
780 ELSE IF (row_indices(irow) > col_indices(icol)) THEN
781 icol = icol + 1
782 ELSE ! diagonal element
783 determinant = determinant*a(irow, icol)
784 irow = irow + 1
785 icol = icol + 1
786 END IF
787 END DO
788 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 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 DEALLOCATE (ipivot)
812 CALL timestop(handle)
813
814 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 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 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
829 COMPLEX(kind=dp), DIMENSION(1) :: work1
830 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: mat
831 INTEGER :: handle, info, lwork, nrows_global
832 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
833
834#if defined(__parallel)
835 INTEGER :: liwork
836 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
837 INTEGER, DIMENSION(1) :: iwork1
838 INTEGER, DIMENSION(9) :: desca
839#else
840 INTEGER :: lda
841#endif
842
843 CALL timeset(routinen, handle)
844
845 mat => matrix%local_data
846 nrows_global = matrix%matrix_struct%nrow_global
847 cpassert(nrows_global == matrix%matrix_struct%ncol_global)
848 ALLOCATE (ipivot(nrows_global))
849
850 ! do LU decomposition
851#if defined(__parallel)
852 desca = matrix%matrix_struct%descriptor
853 CALL pzgetrf(nrows_global, nrows_global, &
854 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 IF (info /= 0) THEN
861 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 ipivot, work1, -1, iwork1, -1, info)
868 lwork = int(work1(1))
869 liwork = int(iwork1(1))
870 ALLOCATE (work(lwork))
871 ALLOCATE (iwork(liwork))
872 CALL pzgetri(nrows_global, mat(1, 1), 1, 1, desca, &
873 ipivot, work, lwork, iwork, liwork, info)
874 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 DEALLOCATE (work)
882 DEALLOCATE (ipivot)
883
884 IF (PRESENT(info_out)) THEN
885 info_out = info
886 ELSE
887 IF (info /= 0) &
888 CALL cp_abort(__location__, "LU inversion has failed")
889 END IF
890
891 CALL timestop(handle)
892
893 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 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 CALL timeset(routinen, handle)
919
920 context => matrix_a%matrix_struct%context
921 myprow = context%mepos(1)
922 mypcol = context%mepos(2)
923 nprow = context%num_pe(1)
924 npcol = context%num_pe(2)
925
926 group = matrix_a%matrix_struct%para_env
927
928 nrow_local = min(matrix_a%matrix_struct%nrow_locals(myprow), matrix_b%matrix_struct%nrow_locals(myprow))
929 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 matrix_b%local_data(1:nrow_local, 1:ncol_local))
934
935 CALL group%sum(trace)
936
937 CALL timestop(handle)
938
939 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 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 CALL timeset(routinen, handle)
994 side_char = 'L'
995 unit_diag = 'N'
996 uplo = 'U'
997 transa = 'N'
998 invert = .false.
999 al = z_one
1000 CALL cp_cfm_get_info(matrix_b, nrow_global=m, ncol_global=n)
1001 IF (PRESENT(side)) side_char = side
1002 IF (PRESENT(invert_tr)) invert = invert_tr
1003 IF (PRESENT(uplo_tr)) uplo = uplo_tr
1004 IF (PRESENT(unit_diag_tr)) THEN
1005 IF (unit_diag_tr) THEN
1006 unit_diag = 'U'
1007 ELSE
1008 unit_diag = 'N'
1009 END IF
1010 END IF
1011 IF (PRESENT(transa_tr)) transa = transa_tr
1012 IF (PRESENT(alpha)) al = alpha
1013 IF (PRESENT(n_rows)) m = n_rows
1014 IF (PRESENT(n_cols)) n = n_cols
1015
1016 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 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 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 CALL timestop(handle)
1049
1050 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 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 POINTER :: a
1070#if defined(__parallel)
1071 INTEGER, DIMENSION(9) :: desca
1072#endif
1073
1074 CALL timeset(routinen, handle)
1075
1076 unit_diag = 'N'
1077 my_uplo = 'U'
1078 IF (PRESENT(uplo)) my_uplo = uplo
1079
1080 ncol_global = matrix_a%matrix_struct%ncol_global
1081
1082 a => matrix_a%local_data
1083
1084#if defined(__parallel)
1085 desca(:) = matrix_a%matrix_struct%descriptor(:)
1086 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 IF (PRESENT(info_out)) THEN
1092 info_out = info
1093 ELSE
1094 IF (info /= 0) &
1095 CALL cp_abort(__location__, &
1096 "triangular invert failed: matrix is not positive definite or ill-conditioned")
1097 END IF
1098
1099 CALL timestop(handle)
1100 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 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 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 CALL timeset(routinen, handle)
1125
1126 nrow_global = matrix%matrix_struct%nrow_global
1127 ncol_global = matrix%matrix_struct%ncol_global
1128
1129 cpassert(matrixt%matrix_struct%nrow_global == ncol_global)
1130 cpassert(matrixt%matrix_struct%ncol_global == nrow_global)
1131
1132 aa => matrix%local_data
1133 cc => matrixt%local_data
1134
1135#if defined(__parallel)
1136 desca = matrix%matrix_struct%descriptor
1137 descc = matrixt%matrix_struct%descriptor
1138 SELECT CASE (trans)
1139 CASE ('T')
1140 CALL pztranu(nrow_global, ncol_global, &
1141 z_one, aa(1, 1), 1, 1, desca, &
1142 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 z_zero, cc(1, 1), 1, 1, descc)
1147 CASE DEFAULT
1148 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 CALL timestop(handle)
1172 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 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 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa
1192 INTEGER :: handle, lwork, ncols, ncols_local, &
1193 nrows, nrows_local
1194 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 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 ncol_local=ncols_local)
1209 aa => matrix%local_data
1210
1211 SELECT CASE (mode)
1212 CASE ('M', 'm')
1213 lwork = 1
1214 CASE ('1', 'O', 'o')
1215#if defined(__parallel)
1216 lwork = ncols_local
1217#else
1218 lwork = 1
1219#endif
1220 CASE ('I', 'i')
1221#if defined(__parallel)
1222 lwork = nrows_local
1223#else
1224 lwork = nrows
1225#endif
1226 CASE ('F', 'f', 'E', 'e')
1227 lwork = 1
1228 CASE DEFAULT
1229 cpabort("mode input is not valid")
1230 END SELECT
1231
1232 ALLOCATE (work(lwork))
1233
1234#if defined(__parallel)
1235 desca = matrix%matrix_struct%descriptor
1236 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 DEALLOCATE (work)
1243 CALL timestop(handle)
1244 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 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 REAL(dp), DIMENSION(:), ALLOCATABLE :: work
1268#endif
1269 CALL timeset(routinen, handle)
1270 CALL cp_cfm_get_info(matrix, ncol_global=ncol)
1271 sn_cmplx = cmplx(sn, 0.0_dp, dp)
1272#if defined(__parallel)
1273 IF (1 /= matrix%matrix_struct%context%n_pid) THEN
1274 lwork = 2*ncol + 1
1275 ALLOCATE (work(lwork))
1276 desc(:) = matrix%matrix_struct%descriptor(:)
1277 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 cs, sn_cmplx, work, lwork, info)
1282 cpassert(info == 0)
1283 DEALLOCATE (work)
1284 ELSE
1285#endif
1286 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 CALL timestop(handle)
1291 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 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 REAL(dp), DIMENSION(:), ALLOCATABLE :: work
1315#endif
1316 CALL timeset(routinen, handle)
1317 CALL cp_cfm_get_info(matrix, nrow_global=nrow)
1318 sn_cmplx = cmplx(sn, 0.0_dp, dp)
1319#if defined(__parallel)
1320 IF (1 /= matrix%matrix_struct%context%n_pid) THEN
1321 lwork = 2*nrow + 1
1322 ALLOCATE (work(lwork))
1323 desc(:) = matrix%matrix_struct%descriptor(:)
1324 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 cs, sn_cmplx, work, lwork, info)
1329 cpassert(info == 0)
1330 DEALLOCATE (work)
1331 ELSE
1332#endif
1333 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 CALL timestop(handle)
1338 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1361 TYPE(cp_cfm_type) :: work
1362
1363 CALL timeset(routinen, handle)
1364
1365 IF (.NOT. PRESENT(workspace)) THEN
1366 CALL cp_cfm_create(work, matrix%matrix_struct)
1367 ELSE
1368 work = workspace
1369 END IF
1370
1371 myuplo = 'U'
1372 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 col_indices=col_indices)
1380
1381 DO jjb = 1, ncol_local
1382 j_global = col_indices(jjb)
1383 DO iib = 1, nrow_local
1384 i_global = row_indices(iib)
1385 IF (merge(j_global < i_global, j_global > i_global, (myuplo == "U") .OR. (myuplo == "u"))) THEN
1386 matrix%local_data(iib, jjb) = z_zero
1387 ELSE IF (j_global == i_global) THEN
1388 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 CALL cp_cfm_transpose(matrix, 'C', work)
1394
1395 CALL cp_cfm_scale_and_add(z_one, matrix, z_one, work)
1396
1397 IF (.NOT. PRESENT(workspace)) THEN
1398 CALL cp_cfm_release(work)
1399 END IF
1400
1401 CALL timestop(handle)
1402
1403 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices
1419
1420 CALL timeset(routinen, handle)
1421
1422 CALL cp_cfm_get_info(matrix, col_indices=col_indices, nrow_local=nrow_local, &
1423 ncol_local=ncol_local)
1424
1425 ! the efficiency could be improved by making use of the row-col distribution of scalapack
1426 norm_array = 0.0_dp
1427 DO j = 1, ncol_local
1428 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 aimag(matrix%local_data(i, j))**2
1432 END DO
1433 END DO
1434 CALL matrix%matrix_struct%para_env%sum(norm_array)
1435 norm_array = sqrt(norm_array)
1436
1437 CALL timestop(handle)
1438 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 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_row, j_col, max_index, ncol_local, &
1452 nrow_local
1453 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1454
1455 max_index = huge(max_index)
1456 IF (PRESENT(n_active)) max_index = n_active
1457
1458 CALL cp_cfm_get_info(matrix=matrix, nrow_local=nrow_local, ncol_local=ncol_local, &
1459 row_indices=row_indices, col_indices=col_indices)
1460
1461 DO j_col = 1, ncol_local
1462 DO i_row = 1, nrow_local
1463 IF (row_indices(i_row) == col_indices(j_col) .AND. row_indices(i_row) <= max_index) THEN
1464 matrix%local_data(i_row, j_col) = matrix%local_data(i_row, j_col) + alpha
1465 END IF
1466 END DO
1467 END DO
1468
1469 END SUBROUTINE cp_cfm_add_on_diag
1470
1471! **************************************************************************************************
1472!> \brief returns the diagonal of a complex full matrix: diag(i)= A_{ii}.
1473!> Each diagonal entry is owned by one process. The sum over the
1474!> process grid collects the entries.
1475!> \param matrix ...
1476!> \param diag ...
1477! **************************************************************************************************
1478 SUBROUTINE cp_cfm_get_diag(matrix, diag)
1479 TYPE(cp_cfm_type), INTENT(IN) :: matrix
1480 COMPLEX(KIND=dp), DIMENSION(:), INTENT(OUT) :: diag
1481
1482 CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_cfm_get_diag'
1483
1484 INTEGER :: handle, i, j, ncol_local, nrow_local
1485 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1486
1487 CALL timeset(routinen, handle)
1488
1489 CALL cp_cfm_get_info(matrix, col_indices=col_indices, row_indices=row_indices, &
1490 nrow_local=nrow_local, ncol_local=ncol_local)
1491
1492 diag = z_zero
1493 DO j = 1, ncol_local
1494 DO i = 1, nrow_local
1495 IF (row_indices(i) == col_indices(j)) THEN
1496 diag(col_indices(j)) = matrix%local_data(i, j)
1497 END IF
1498 END DO
1499 END DO
1500 CALL matrix%matrix_struct%para_env%sum(diag)
1501
1502 CALL timestop(handle)
1503 END SUBROUTINE cp_cfm_get_diag
1504
1505END MODULE cp_cfm_basic_linalg
cp_cfm_basic_linalg::cp_cfm_scale
Definition cp_cfm_basic_linalg.F:63
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:57
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
kahan_sum::accurate_dot_product
Definition kahan_sum.F:52
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_scale_and_add
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
Definition cp_cfm_basic_linalg.F:243
cp_cfm_basic_linalg::cp_cfm_lu_invert
subroutine, public cp_cfm_lu_invert(matrix, info_out)
Inverts a matrix using LU decomposition. The input matrix will be overwritten.
Definition cp_cfm_basic_linalg.F:823
cp_cfm_basic_linalg::cp_cfm_get_diag
subroutine, public cp_cfm_get_diag(matrix, diag)
returns the diagonal of a complex full matrix: diag(i)= A_{ii}. Each diagonal entry is owned by one p...
Definition cp_cfm_basic_linalg.F:1479
cp_cfm_basic_linalg::cp_cfm_norm
real(kind=dp) function, public cp_cfm_norm(matrix, mode)
Norm of matrix using (p)zlange.
Definition cp_cfm_basic_linalg.F:1185
cp_cfm_basic_linalg::cp_cfm_gemm
subroutine, public cp_cfm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, c_first_row)
Performs one of the matrix-matrix operations: matrix_c = alpha * op1( matrix_a ) * op2( matrix_b ) + ...
Definition cp_cfm_basic_linalg.F:562
cp_cfm_basic_linalg::cp_cfm_solve
subroutine, public cp_cfm_solve(matrix_a, general_a, determinant)
Solve the system of linear equations A*b=A_general using LU decomposition. Pay attention that both ma...
Definition cp_cfm_basic_linalg.F:731
cp_cfm_basic_linalg::cp_cfm_transpose
subroutine, public cp_cfm_transpose(matrix, trans, matrixt)
Transposes a BLACS distributed complex matrix.
Definition cp_cfm_basic_linalg.F:1110
cp_cfm_basic_linalg::cp_cfm_rot_rows
subroutine, public cp_cfm_rot_rows(matrix, irow, jrow, cs, sn)
Applies a planar rotation defined by cs and sn to the i'th and j'th rows.
Definition cp_cfm_basic_linalg.F:1256
cp_cfm_basic_linalg::cp_cfm_scale_and_add_fm
subroutine, public cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b). where b is a real matrix (adapted from cp_cf...
Definition cp_cfm_basic_linalg.F:352
cp_cfm_basic_linalg::cp_cfm_schur_product
subroutine, public cp_cfm_schur_product(matrix_a, matrix_b, matrix_c)
Computes the element-wise (Schur) product of two matrices: C = A \circ B .
Definition cp_cfm_basic_linalg.F:156
cp_cfm_basic_linalg::cp_cfm_triangular_multiply
subroutine, public cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
Multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
Definition cp_cfm_basic_linalg.F:978
cp_cfm_basic_linalg::cp_cfm_vectorsnorm
subroutine, public cp_cfm_vectorsnorm(matrix, norm_array)
find the norm of each column norm_{j}= sqrt( \sum_{i} A_{ij}*conjg(A_{ij}) ) Complex-valued mirror of...
Definition cp_cfm_basic_linalg.F:1412
cp_cfm_basic_linalg::cp_cfm_uplo_to_full
subroutine, public cp_cfm_uplo_to_full(matrix, workspace, uplo)
...
Definition cp_cfm_basic_linalg.F:1350
cp_cfm_basic_linalg::cp_cfm_add_on_diag
subroutine, public cp_cfm_add_on_diag(matrix, alpha, n_active)
Adds a scalar to the diagonal of a distributed complex full matrix.
Definition cp_cfm_basic_linalg.F:1447
cp_cfm_basic_linalg::cp_cfm_det
subroutine, public cp_cfm_det(matrix_a, det_a)
Computes the determinant (with a correct sign even in parallel environment!) of a complex square matr...
Definition cp_cfm_basic_linalg.F:78
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling)
Scales columns of the full matrix by corresponding factors.
Definition cp_cfm_basic_linalg.F:632
cp_cfm_basic_linalg::cp_cfm_rot_cols
subroutine, public cp_cfm_rot_cols(matrix, icol, jcol, cs, sn)
Applies a planar rotation defined by cs and sn to the i'th and j'th columnns.
Definition cp_cfm_basic_linalg.F:1303
cp_cfm_basic_linalg::cp_cfm_triangular_invert
subroutine, public cp_cfm_triangular_invert(matrix_a, uplo, info_out)
Inverts a triangular matrix.
Definition cp_cfm_basic_linalg.F:1060
cp_cfm_basic_linalg::cp_cfm_lu_decompose
subroutine, public cp_cfm_lu_decompose(matrix_a, determinant)
Computes LU decomposition of a given matrix.
Definition cp_cfm_basic_linalg.F:457
cp_cfm_basic_linalg::cp_cfm_trace
subroutine, public cp_cfm_trace(matrix_a, matrix_b, trace)
Returns the trace of matrix_a^T matrix_b, i.e sum_{i,j}(matrix_a(i,j)*matrix_b(i,j)) .
Definition cp_cfm_basic_linalg.F:908
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:183
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:126
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:632
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_equivalent
logical function, public cp_fm_struct_equivalent(fmstruct1, fmstruct2)
returns true if the two matrix structures are equivalent, false otherwise.
Definition cp_fm_struct.F:392
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:144
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:70
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
message_passing::mp_comm_type
Definition message_passing.F:179