(git:858b7a1)
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, size_a
710
711 CALL timeset(routinen, handle)
712
713 NULLIFY (a)
714
715 a => matrix_a%local_data
716 size_a = SIZE(a, 1)*SIZE(a, 2)
717
718 CALL zscal(size_a, alpha, a(1, 1), 1)
719
720 CALL timestop(handle)
721 END SUBROUTINE cp_cfm_zscale
722
723! **************************************************************************************************
724!> \brief Solve the system of linear equations A*b=A_general using LU decomposition.
725!> Pay attention that both matrices are overwritten on exit and that
726!> the result is stored into the matrix 'general_a'.
727!> \param matrix_a matrix A (overwritten on exit)
728!> \param general_a (input) matrix A_general, (output) matrix B
729!> \param determinant (optional) determinant
730!> \author Florian Schiffmann
731! **************************************************************************************************
732 SUBROUTINE cp_cfm_solve(matrix_a, general_a, determinant)
733 TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, general_a
734 COMPLEX(kind=dp), OPTIONAL :: determinant
735
736 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_solve'
737
738 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, a_general
739 INTEGER :: counter, handle, info, irow, nrow_global
740 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
741
742#if defined(__parallel)
743 INTEGER :: icol, ncol_local, nrow_local
744 INTEGER, DIMENSION(9) :: desca, descb
745 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
746#else
747 INTEGER :: lda, ldb
748#endif
749
750 CALL timeset(routinen, handle)
751
752 a => matrix_a%local_data
753 a_general => general_a%local_data
754 nrow_global = matrix_a%matrix_struct%nrow_global
755 ALLOCATE (ipivot(nrow_global))
756
757#if defined(__parallel)
758 desca(:) = matrix_a%matrix_struct%descriptor(:)
759 descb(:) = general_a%matrix_struct%descriptor(:)
760 CALL pzgetrf(nrow_global, nrow_global, a(1, 1), 1, 1, desca, ipivot, info)
761 IF (PRESENT(determinant)) THEN
762 CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, &
763 row_indices=row_indices, col_indices=col_indices)
764
765 counter = 0
766 DO irow = 1, nrow_local
767 IF (ipivot(irow) .NE. row_indices(irow)) counter = counter + 1
768 END DO
769
770 IF (mod(counter, 2) == 0) THEN
771 determinant = z_one
772 ELSE
773 determinant = -z_one
774 END IF
775
776 ! compute product of diagonal elements
777 irow = 1
778 icol = 1
779 DO WHILE (irow <= nrow_local .AND. icol <= ncol_local)
780 IF (row_indices(irow) < col_indices(icol)) THEN
781 irow = irow + 1
782 ELSE IF (row_indices(irow) > col_indices(icol)) THEN
783 icol = icol + 1
784 ELSE ! diagonal element
785 determinant = determinant*a(irow, icol)
786 irow = irow + 1
787 icol = icol + 1
788 END IF
789 END DO
790 CALL matrix_a%matrix_struct%para_env%prod(determinant)
791 END IF
792
793 CALL pzgetrs("N", nrow_global, nrow_global, a(1, 1), 1, 1, desca, &
794 ipivot, a_general(1, 1), 1, 1, descb, info)
795#else
796 lda = SIZE(a, 1)
797 ldb = SIZE(a_general, 1)
798 CALL zgetrf(nrow_global, nrow_global, a(1, 1), lda, ipivot, info)
799 IF (PRESENT(determinant)) THEN
800 counter = 0
801 determinant = z_one
802 DO irow = 1, nrow_global
803 IF (ipivot(irow) .NE. irow) counter = counter + 1
804 determinant = determinant*a(irow, irow)
805 END DO
806 IF (mod(counter, 2) == 1) determinant = -1.0_dp*determinant
807 END IF
808 CALL zgetrs("N", nrow_global, nrow_global, a(1, 1), lda, ipivot, a_general(1, 1), ldb, info)
809#endif
810
811 ! info is allowed to be zero
812 ! this does just signal a zero diagonal element
813 DEALLOCATE (ipivot)
814 CALL timestop(handle)
815
816 END SUBROUTINE cp_cfm_solve
817
818! **************************************************************************************************
819!> \brief Inverts a matrix using LU decomposition. The input matrix will be overwritten.
820!> \param matrix input a general square non-singular matrix, outputs its inverse
821!> \param info_out optional, if present outputs the info from (p)zgetri
822!> \author Lianheng Tong
823! **************************************************************************************************
824 SUBROUTINE cp_cfm_lu_invert(matrix, info_out)
825 TYPE(cp_cfm_type), INTENT(IN) :: matrix
826 INTEGER, INTENT(out), OPTIONAL :: info_out
827
828 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_lu_invert'
829
830 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
831 COMPLEX(kind=dp), DIMENSION(1) :: work1
832 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: mat
833 INTEGER :: handle, info, lwork, nrows_global
834 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
835
836#if defined(__parallel)
837 INTEGER :: liwork
838 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
839 INTEGER, DIMENSION(1) :: iwork1
840 INTEGER, DIMENSION(9) :: desca
841#else
842 INTEGER :: lda
843#endif
844
845 CALL timeset(routinen, handle)
846
847 mat => matrix%local_data
848 nrows_global = matrix%matrix_struct%nrow_global
849 cpassert(nrows_global .EQ. matrix%matrix_struct%ncol_global)
850 ALLOCATE (ipivot(nrows_global))
851
852 ! do LU decomposition
853#if defined(__parallel)
854 desca = matrix%matrix_struct%descriptor
855 CALL pzgetrf(nrows_global, nrows_global, &
856 mat(1, 1), 1, 1, desca, ipivot, info)
857#else
858 lda = SIZE(mat, 1)
859 CALL zgetrf(nrows_global, nrows_global, &
860 mat(1, 1), lda, ipivot, info)
861#endif
862 IF (info /= 0) THEN
863 CALL cp_abort(__location__, "LU decomposition has failed")
864 END IF
865
866 ! do inversion
867#if defined(__parallel)
868 CALL pzgetri(nrows_global, mat(1, 1), 1, 1, desca, &
869 ipivot, work1, -1, iwork1, -1, info)
870 lwork = int(work1(1))
871 liwork = int(iwork1(1))
872 ALLOCATE (work(lwork))
873 ALLOCATE (iwork(liwork))
874 CALL pzgetri(nrows_global, mat(1, 1), 1, 1, desca, &
875 ipivot, work, lwork, iwork, liwork, info)
876 DEALLOCATE (iwork)
877#else
878 CALL zgetri(nrows_global, mat(1, 1), lda, ipivot, work1, -1, info)
879 lwork = int(work1(1))
880 ALLOCATE (work(lwork))
881 CALL zgetri(nrows_global, mat(1, 1), lda, ipivot, work, lwork, info)
882#endif
883 DEALLOCATE (work)
884 DEALLOCATE (ipivot)
885
886 IF (PRESENT(info_out)) THEN
887 info_out = info
888 ELSE
889 IF (info /= 0) &
890 CALL cp_abort(__location__, "LU inversion has failed")
891 END IF
892
893 CALL timestop(handle)
894
895 END SUBROUTINE cp_cfm_lu_invert
896
897! **************************************************************************************************
898!> \brief Returns the trace of matrix_a^T matrix_b, i.e
899!> sum_{i,j}(matrix_a(i,j)*matrix_b(i,j)) .
900!> \param matrix_a a complex matrix
901!> \param matrix_b another complex matrix
902!> \param trace value of the trace operator
903!> \par History
904!> * 09.2017 created [Sergey Chulkov]
905!> \author Sergey Chulkov
906!> \note
907!> Based on the subroutine cp_fm_trace(). Note the transposition of matrix_a!
908! **************************************************************************************************
909 SUBROUTINE cp_cfm_trace(matrix_a, matrix_b, trace)
910 TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, matrix_b
911 COMPLEX(kind=dp), INTENT(out) :: trace
912
913 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_trace'
914
915 INTEGER :: handle, mypcol, myprow, ncol_local, &
916 npcol, nprow, nrow_local
917 TYPE(cp_blacs_env_type), POINTER :: context
918 TYPE(mp_comm_type) :: group
919
920 CALL timeset(routinen, handle)
921
922 context => matrix_a%matrix_struct%context
923 myprow = context%mepos(1)
924 mypcol = context%mepos(2)
925 nprow = context%num_pe(1)
926 npcol = context%num_pe(2)
927
928 group = matrix_a%matrix_struct%para_env
929
930 nrow_local = min(matrix_a%matrix_struct%nrow_locals(myprow), matrix_b%matrix_struct%nrow_locals(myprow))
931 ncol_local = min(matrix_a%matrix_struct%ncol_locals(mypcol), matrix_b%matrix_struct%ncol_locals(mypcol))
932
933 ! compute an accurate dot-product
934 trace = accurate_dot_product(matrix_a%local_data(1:nrow_local, 1:ncol_local), &
935 matrix_b%local_data(1:nrow_local, 1:ncol_local))
936
937 CALL group%sum(trace)
938
939 CALL timestop(handle)
940
941 END SUBROUTINE cp_cfm_trace
942
943! **************************************************************************************************
944!> \brief Multiplies in place by a triangular matrix:
945!> matrix_b = alpha op(triangular_matrix) matrix_b
946!> or (if side='R')
947!> matrix_b = alpha matrix_b op(triangular_matrix)
948!> op(triangular_matrix) is:
949!> triangular_matrix (if transa="N" and invert_tr=.false.)
950!> triangular_matrix^T (if transa="T" and invert_tr=.false.)
951!> triangular_matrix^H (if transa="C" and invert_tr=.false.)
952!> triangular_matrix^(-1) (if transa="N" and invert_tr=.true.)
953!> triangular_matrix^(-T) (if transa="T" and invert_tr=.true.)
954!> triangular_matrix^(-H) (if transa="C" and invert_tr=.true.)
955!> \param triangular_matrix the triangular matrix that multiplies the other
956!> \param matrix_b the matrix that gets multiplied and stores the result
957!> \param side on which side of matrix_b stays op(triangular_matrix)
958!> (defaults to 'L')
959!> \param transa_tr ...
960!> \param invert_tr if the triangular matrix should be inverted
961!> (defaults to false)
962!> \param uplo_tr if triangular_matrix is stored in the upper ('U') or
963!> lower ('L') triangle (defaults to 'U')
964!> \param unit_diag_tr if the diagonal elements of triangular_matrix should
965!> be assumed to be 1 (defaults to false)
966!> \param n_rows the number of rows of the result (defaults to
967!> size(matrix_b,1))
968!> \param n_cols the number of columns of the result (defaults to
969!> size(matrix_b,2))
970!> \param alpha ...
971!> \par History
972!> 08.2002 created [fawzi]
973!> \author Fawzi Mohamed
974!> \note
975!> needs an mpi env
976! **************************************************************************************************
977 SUBROUTINE cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, &
978 transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, &
979 alpha)
980 TYPE(cp_cfm_type), INTENT(IN) :: triangular_matrix, matrix_b
981 CHARACTER, INTENT(in), OPTIONAL :: side, transa_tr
982 LOGICAL, INTENT(in), OPTIONAL :: invert_tr
983 CHARACTER, INTENT(in), OPTIONAL :: uplo_tr
984 LOGICAL, INTENT(in), OPTIONAL :: unit_diag_tr
985 INTEGER, INTENT(in), OPTIONAL :: n_rows, n_cols
986 COMPLEX(kind=dp), INTENT(in), OPTIONAL :: alpha
987
988 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_triangular_multiply'
989
990 CHARACTER :: side_char, transa, unit_diag, uplo
991 COMPLEX(kind=dp) :: al
992 INTEGER :: handle, m, n
993 LOGICAL :: invert
994
995 CALL timeset(routinen, handle)
996 side_char = 'L'
997 unit_diag = 'N'
998 uplo = 'U'
999 transa = 'N'
1000 invert = .false.
1001 al = cmplx(1.0_dp, 0.0_dp, dp)
1002 CALL cp_cfm_get_info(matrix_b, nrow_global=m, ncol_global=n)
1003 IF (PRESENT(side)) side_char = side
1004 IF (PRESENT(invert_tr)) invert = invert_tr
1005 IF (PRESENT(uplo_tr)) uplo = uplo_tr
1006 IF (PRESENT(unit_diag_tr)) THEN
1007 IF (unit_diag_tr) THEN
1008 unit_diag = 'U'
1009 ELSE
1010 unit_diag = 'N'
1011 END IF
1012 END IF
1013 IF (PRESENT(transa_tr)) transa = transa_tr
1014 IF (PRESENT(alpha)) al = alpha
1015 IF (PRESENT(n_rows)) m = n_rows
1016 IF (PRESENT(n_cols)) n = n_cols
1017
1018 IF (invert) THEN
1019
1020#if defined(__parallel)
1021 CALL pztrsm(side_char, uplo, transa, unit_diag, m, n, al, &
1022 triangular_matrix%local_data(1, 1), 1, 1, &
1023 triangular_matrix%matrix_struct%descriptor, &
1024 matrix_b%local_data(1, 1), 1, 1, &
1025 matrix_b%matrix_struct%descriptor(1))
1026#else
1027 CALL ztrsm(side_char, uplo, transa, unit_diag, m, n, al, &
1028 triangular_matrix%local_data(1, 1), &
1029 SIZE(triangular_matrix%local_data, 1), &
1030 matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1))
1031#endif
1032
1033 ELSE
1034
1035#if defined(__parallel)
1036 CALL pztrmm(side_char, uplo, transa, unit_diag, m, n, al, &
1037 triangular_matrix%local_data(1, 1), 1, 1, &
1038 triangular_matrix%matrix_struct%descriptor, &
1039 matrix_b%local_data(1, 1), 1, 1, &
1040 matrix_b%matrix_struct%descriptor(1))
1041#else
1042 CALL ztrmm(side_char, uplo, transa, unit_diag, m, n, al, &
1043 triangular_matrix%local_data(1, 1), &
1044 SIZE(triangular_matrix%local_data, 1), &
1045 matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1))
1046#endif
1047
1048 END IF
1049
1050 CALL timestop(handle)
1051
1052 END SUBROUTINE cp_cfm_triangular_multiply
1053
1054! **************************************************************************************************
1055!> \brief Inverts a triangular matrix.
1056!> \param matrix_a ...
1057!> \param uplo ...
1058!> \param info_out ...
1059!> \author MI
1060! **************************************************************************************************
1061 SUBROUTINE cp_cfm_triangular_invert(matrix_a, uplo, info_out)
1062 TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
1063 CHARACTER, INTENT(in), OPTIONAL :: uplo
1064 INTEGER, INTENT(out), OPTIONAL :: info_out
1065
1066 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_triangular_invert'
1067
1068 CHARACTER :: unit_diag, my_uplo
1069 INTEGER :: handle, info, ncol_global
1070 COMPLEX(kind=dp), DIMENSION(:, :), &
1071 POINTER :: a
1072#if defined(__parallel)
1073 INTEGER, DIMENSION(9) :: desca
1074#endif
1075
1076 CALL timeset(routinen, handle)
1077
1078 unit_diag = 'N'
1079 my_uplo = 'U'
1080 IF (PRESENT(uplo)) my_uplo = uplo
1081
1082 ncol_global = matrix_a%matrix_struct%ncol_global
1083
1084 a => matrix_a%local_data
1085
1086#if defined(__parallel)
1087 desca(:) = matrix_a%matrix_struct%descriptor(:)
1088 CALL pztrtri(my_uplo, unit_diag, ncol_global, a(1, 1), 1, 1, desca, info)
1089#else
1090 CALL ztrtri(my_uplo, unit_diag, ncol_global, a(1, 1), ncol_global, info)
1091#endif
1092
1093 IF (PRESENT(info_out)) THEN
1094 info_out = info
1095 ELSE
1096 IF (info /= 0) &
1097 CALL cp_abort(__location__, &
1098 "triangular invert failed: matrix is not positive definite or ill-conditioned")
1099 END IF
1100
1101 CALL timestop(handle)
1102 END SUBROUTINE cp_cfm_triangular_invert
1103
1104! **************************************************************************************************
1105!> \brief Transposes a BLACS distributed complex matrix.
1106!> \param matrix input matrix
1107!> \param trans 'T' for transpose, 'C' for Hermitian conjugate
1108!> \param matrixt output matrix
1109!> \author Lianheng Tong
1110! **************************************************************************************************
1111 SUBROUTINE cp_cfm_transpose(matrix, trans, matrixt)
1112 TYPE(cp_cfm_type), INTENT(IN) :: matrix
1113 CHARACTER, INTENT(in) :: trans
1114 TYPE(cp_cfm_type), INTENT(IN) :: matrixt
1115
1116 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_transpose'
1117
1118 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa, cc
1119 INTEGER :: handle, ncol_global, nrow_global
1120#if defined(__parallel)
1121 INTEGER, DIMENSION(9) :: desca, descc
1122#elif !defined(__MKL)
1123 INTEGER :: ii, jj
1124#endif
1125
1126 CALL timeset(routinen, handle)
1127
1128 nrow_global = matrix%matrix_struct%nrow_global
1129 ncol_global = matrix%matrix_struct%ncol_global
1130
1131 cpassert(matrixt%matrix_struct%nrow_global == ncol_global)
1132 cpassert(matrixt%matrix_struct%ncol_global == nrow_global)
1133
1134 aa => matrix%local_data
1135 cc => matrixt%local_data
1136
1137#if defined(__parallel)
1138 desca = matrix%matrix_struct%descriptor
1139 descc = matrixt%matrix_struct%descriptor
1140 SELECT CASE (trans)
1141 CASE ('T')
1142 CALL pztranu(nrow_global, ncol_global, &
1143 z_one, aa(1, 1), 1, 1, desca, &
1144 z_zero, cc(1, 1), 1, 1, descc)
1145 CASE ('C')
1146 CALL pztranc(nrow_global, ncol_global, &
1147 z_one, aa(1, 1), 1, 1, desca, &
1148 z_zero, cc(1, 1), 1, 1, descc)
1149 CASE DEFAULT
1150 cpabort("trans only accepts 'T' or 'C'")
1151 END SELECT
1152#elif defined(__MKL)
1153 CALL mkl_zomatcopy('C', trans, nrow_global, ncol_global, 1.0_dp, aa(1, 1), nrow_global, cc(1, 1), ncol_global)
1154#else
1155 SELECT CASE (trans)
1156 CASE ('T')
1157 DO jj = 1, ncol_global
1158 DO ii = 1, nrow_global
1159 cc(ii, jj) = aa(jj, ii)
1160 END DO
1161 END DO
1162 CASE ('C')
1163 DO jj = 1, ncol_global
1164 DO ii = 1, nrow_global
1165 cc(ii, jj) = conjg(aa(jj, ii))
1166 END DO
1167 END DO
1168 CASE DEFAULT
1169 cpabort("trans only accepts 'T' or 'C'")
1170 END SELECT
1171#endif
1172
1173 CALL timestop(handle)
1174 END SUBROUTINE cp_cfm_transpose
1175
1176! **************************************************************************************************
1177!> \brief Norm of matrix using (p)zlange.
1178!> \param matrix input a general matrix
1179!> \param mode 'M' max abs element value,
1180!> '1' or 'O' one norm, i.e. maximum column sum,
1181!> 'I' infinity norm, i.e. maximum row sum,
1182!> 'F' or 'E' Frobenius norm, i.e. sqrt of sum of all squares of elements
1183!> \return the norm according to mode
1184!> \author Lianheng Tong
1185! **************************************************************************************************
1186 FUNCTION cp_cfm_norm(matrix, mode) RESULT(res)
1187 TYPE(cp_cfm_type), INTENT(IN) :: matrix
1188 CHARACTER, INTENT(IN) :: mode
1189 REAL(kind=dp) :: res
1190
1191 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_norm'
1192
1193 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa
1194 INTEGER :: handle, lwork, ncols, ncols_local, &
1195 nrows, nrows_local
1196 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
1197
1198#if defined(__parallel)
1199 INTEGER, DIMENSION(9) :: desca
1200#else
1201 INTEGER :: lda
1202#endif
1203
1204 CALL timeset(routinen, handle)
1205
1206 CALL cp_cfm_get_info(matrix=matrix, &
1207 nrow_global=nrows, &
1208 ncol_global=ncols, &
1209 nrow_local=nrows_local, &
1210 ncol_local=ncols_local)
1211 aa => matrix%local_data
1212
1213 SELECT CASE (mode)
1214 CASE ('M', 'm')
1215 lwork = 1
1216 CASE ('1', 'O', 'o')
1217#if defined(__parallel)
1218 lwork = ncols_local
1219#else
1220 lwork = 1
1221#endif
1222 CASE ('I', 'i')
1223#if defined(__parallel)
1224 lwork = nrows_local
1225#else
1226 lwork = nrows
1227#endif
1228 CASE ('F', 'f', 'E', 'e')
1229 lwork = 1
1230 CASE DEFAULT
1231 cpabort("mode input is not valid")
1232 END SELECT
1233
1234 ALLOCATE (work(lwork))
1235
1236#if defined(__parallel)
1237 desca = matrix%matrix_struct%descriptor
1238 res = pzlange(mode, nrows, ncols, aa(1, 1), 1, 1, desca, work)
1239#else
1240 lda = SIZE(aa, 1)
1241 res = zlange(mode, nrows, ncols, aa(1, 1), lda, work)
1242#endif
1243
1244 DEALLOCATE (work)
1245 CALL timestop(handle)
1246 END FUNCTION cp_cfm_norm
1247
1248! **************************************************************************************************
1249!> \brief Applies a planar rotation defined by cs and sn to the i'th and j'th rows.
1250!> \param matrix ...
1251!> \param irow ...
1252!> \param jrow ...
1253!> \param cs cosine of the rotation angle
1254!> \param sn sinus of the rotation angle
1255!> \author Ole Schuett
1256! **************************************************************************************************
1257 SUBROUTINE cp_cfm_rot_rows(matrix, irow, jrow, cs, sn)
1258 TYPE(cp_cfm_type), INTENT(IN) :: matrix
1259 INTEGER, INTENT(IN) :: irow, jrow
1260 REAL(dp), INTENT(IN) :: cs, sn
1261
1262 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_rot_rows'
1263 INTEGER :: handle, ncol
1264 COMPLEX(KIND=dp) :: sn_cmplx
1265
1266#if defined(__parallel)
1267 INTEGER :: info, lwork
1268 INTEGER, DIMENSION(9) :: desc
1269 REAL(dp), DIMENSION(:), ALLOCATABLE :: work
1270#endif
1271 CALL timeset(routinen, handle)
1272 CALL cp_cfm_get_info(matrix, ncol_global=ncol)
1273 sn_cmplx = cmplx(sn, 0.0_dp, dp)
1274#if defined(__parallel)
1275 IF (1 /= matrix%matrix_struct%context%n_pid) THEN
1276 lwork = 2*ncol + 1
1277 ALLOCATE (work(lwork))
1278 desc(:) = matrix%matrix_struct%descriptor(:)
1279 CALL pzrot(ncol, &
1280 matrix%local_data(1, 1), irow, 1, desc, ncol, &
1281 matrix%local_data(1, 1), jrow, 1, desc, ncol, &
1282 cs, sn_cmplx, work, lwork, info)
1283 cpassert(info == 0)
1284 DEALLOCATE (work)
1285 ELSE
1286#endif
1287 CALL zrot(ncol, matrix%local_data(irow, 1), ncol, matrix%local_data(jrow, 1), ncol, cs, sn_cmplx)
1288#if defined(__parallel)
1289 END IF
1290#endif
1291 CALL timestop(handle)
1292 END SUBROUTINE cp_cfm_rot_rows
1293
1294! **************************************************************************************************
1295!> \brief Applies a planar rotation defined by cs and sn to the i'th and j'th columnns.
1296!> \param matrix ...
1297!> \param icol ...
1298!> \param jcol ...
1299!> \param cs cosine of the rotation angle
1300!> \param sn sinus of the rotation angle
1301!> \author Ole Schuett
1302! **************************************************************************************************
1303 SUBROUTINE cp_cfm_rot_cols(matrix, icol, jcol, cs, sn)
1304 TYPE(cp_cfm_type), INTENT(IN) :: matrix
1305 INTEGER, INTENT(IN) :: icol, jcol
1306 REAL(dp), INTENT(IN) :: cs, sn
1307
1308 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_rot_cols'
1309 INTEGER :: handle, nrow
1310 COMPLEX(KIND=dp) :: sn_cmplx
1311
1312#if defined(__parallel)
1313 INTEGER :: info, lwork
1314 INTEGER, DIMENSION(9) :: desc
1315 REAL(dp), DIMENSION(:), ALLOCATABLE :: work
1316#endif
1317 CALL timeset(routinen, handle)
1318 CALL cp_cfm_get_info(matrix, nrow_global=nrow)
1319 sn_cmplx = cmplx(sn, 0.0_dp, dp)
1320#if defined(__parallel)
1321 IF (1 /= matrix%matrix_struct%context%n_pid) THEN
1322 lwork = 2*nrow + 1
1323 ALLOCATE (work(lwork))
1324 desc(:) = matrix%matrix_struct%descriptor(:)
1325 CALL pzrot(nrow, &
1326 matrix%local_data(1, 1), 1, icol, desc, 1, &
1327 matrix%local_data(1, 1), 1, jcol, desc, 1, &
1328 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
1405END 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:825
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:1187
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:733
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:1112
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:1258
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:980
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_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:1304
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:1062
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:910
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