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