(git:c835e6b)
Loading...
Searching...
No Matches
cp_cfm_basic_linalg.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Basic linear algebra operations for complex full matrices.
10!> \note
11!> - not all functionality implemented
12!> \par History
13!> Nearly literal copy of Fawzi's routines
14!> \author Joost VandeVondele
15! **************************************************************************************************
16MODULE cp_cfm_basic_linalg
17 USE cp_blacs_env, ONLY: cp_blacs_env_type
18 USE cp_cfm_types, ONLY: cp_cfm_create,&
19 cp_cfm_get_info,&
20 cp_cfm_release,&
21 cp_cfm_to_cfm,&
22 cp_cfm_type
23 USE cp_fm_struct, ONLY: cp_fm_struct_equivalent
24 USE cp_fm_types, ONLY: cp_fm_type
25 USE cp_log_handling, ONLY: cp_to_string
26 USE kahan_sum, ONLY: accurate_dot_product
27 USE kinds, ONLY: dp
28 USE mathconstants, ONLY: z_one,&
29 z_zero
30 USE message_passing, ONLY: mp_comm_type
31#include "../base/base_uses.f90"
32
33 IMPLICIT NONE
34 PRIVATE
35
36 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_basic_linalg'
38
39 PUBLIC :: cp_cfm_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 CALL cp_abort(__location__, &
330 "cp_cfm_scale_and_add is not yet implemented for cases "// &
331 "where input two matrix structures are not equivalent")
332 END IF
333 END IF
334 CALL timestop(handle)
335 END SUBROUTINE cp_cfm_scale_and_add
336
337! **************************************************************************************************
338!> \brief Scale and add two BLACS matrices (a = alpha*a + beta*b).
339!> where b is a real matrix (adapted from cp_cfm_scale_and_add).
340!> \param alpha ...
341!> \param matrix_a ...
342!> \param beta ...
343!> \param matrix_b ...
344!> \date 01.08.2014
345!> \author JGH
346!> \version 1.0
347! **************************************************************************************************
348 SUBROUTINE cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
349 COMPLEX(kind=dp), INTENT(in) :: alpha
350 TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
351 COMPLEX(kind=dp), INTENT(in) :: beta
352 TYPE(cp_fm_type), INTENT(IN) :: matrix_b
353
354 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_scale_and_add_fm'
355
356 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
357 INTEGER :: handle, icol_local, irow_local, mypcol, &
358 myprow, ncol_local, nrow_local
359 REAL(kind=dp), DIMENSION(:, :), POINTER :: b
360
361 CALL timeset(routinen, handle)
362
363 NULLIFY (a, b)
364
365 myprow = matrix_a%matrix_struct%context%mepos(1)
366 mypcol = matrix_a%matrix_struct%context%mepos(2)
367
368 nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
369 ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
370
371 a => matrix_a%local_data
372
373 IF (beta == z_zero) THEN
374
375 IF (alpha == z_zero) THEN
376 a(:, :) = z_zero
377 ELSE IF (alpha == z_one) THEN
378 CALL timestop(handle)
379 RETURN
380 ELSE
381 a(:, :) = alpha*a(:, :)
382 END IF
383
384 ELSE
385 IF (matrix_a%matrix_struct%context /= matrix_b%matrix_struct%context) &
386 cpabort("matrices must be in the same blacs context")
387
388 IF (cp_fm_struct_equivalent(matrix_a%matrix_struct, &
389 matrix_b%matrix_struct)) THEN
390
391 b => matrix_b%local_data
392
393 IF (alpha == z_zero) THEN
394 IF (beta == z_one) THEN
395 !a(:, :) = b(:, :)
396 DO icol_local = 1, ncol_local
397 DO irow_local = 1, nrow_local
398 a(irow_local, icol_local) = b(irow_local, icol_local)
399 END DO
400 END DO
401 ELSE
402 !a(:, :) = beta*b(:, :)
403 DO icol_local = 1, ncol_local
404 DO irow_local = 1, nrow_local
405 a(irow_local, icol_local) = beta*b(irow_local, icol_local)
406 END DO
407 END DO
408 END IF
409 ELSE IF (alpha == z_one) THEN
410 IF (beta == z_one) THEN
411 !a(:, :) = a(:, :)+b(:, :)
412 DO icol_local = 1, ncol_local
413 DO irow_local = 1, nrow_local
414 a(irow_local, icol_local) = a(irow_local, icol_local) + b(irow_local, icol_local)
415 END DO
416 END DO
417 ELSE
418 !a(:, :) = a(:, :)+beta*b(:, :)
419 DO icol_local = 1, ncol_local
420 DO irow_local = 1, nrow_local
421 a(irow_local, icol_local) = a(irow_local, icol_local) + beta*b(irow_local, icol_local)
422 END DO
423 END DO
424 END IF
425 ELSE
426 !a(:, :) = alpha*a(:, :)+beta*b(:, :)
427 DO icol_local = 1, ncol_local
428 DO irow_local = 1, nrow_local
429 a(irow_local, icol_local) = alpha*a(irow_local, icol_local) + beta*b(irow_local, icol_local)
430 END DO
431 END DO
432 END IF
433 ELSE
434 CALL cp_abort(__location__, &
435 "cp_cfm_scale_and_add_fm is not yet implemented for cases "// &
436 "where two input matrix structures are not equivalent")
437 END IF
438 END IF
439 CALL timestop(handle)
440 END SUBROUTINE cp_cfm_scale_and_add_fm
441
442! **************************************************************************************************
443!> \brief Computes LU decomposition of a given matrix.
444!> \param matrix_a full matrix
445!> \param determinant determinant
446!> \date 11.06.2001
447!> \author Matthias Krack
448!> \version 1.0
449!> \note
450!> The actual purpose right now is to efficiently compute the determinant of a given matrix.
451!> The original content of the matrix is destroyed.
452! **************************************************************************************************
453 SUBROUTINE cp_cfm_lu_decompose(matrix_a, determinant)
454 TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
455 COMPLEX(kind=dp), INTENT(out) :: determinant
456
457 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_lu_decompose'
458
459 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
460 INTEGER :: counter, handle, info, irow, nrow_global
461 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
462
463#if defined(__parallel)
464 INTEGER :: icol, ncol_local, nrow_local
465 INTEGER, DIMENSION(9) :: desca
466 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
467#else
468 INTEGER :: lda
469#endif
470
471 CALL timeset(routinen, handle)
472
473 nrow_global = matrix_a%matrix_struct%nrow_global
474 a => matrix_a%local_data
475
476 ALLOCATE (ipivot(nrow_global))
477#if defined(__parallel)
478 CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, &
479 row_indices=row_indices, col_indices=col_indices)
480
481 desca(:) = matrix_a%matrix_struct%descriptor(:)
482 CALL pzgetrf(nrow_global, nrow_global, a(1, 1), 1, 1, desca, ipivot, info)
483
484 counter = 0
485 DO irow = 1, nrow_local
486 IF (ipivot(irow) /= row_indices(irow)) counter = counter + 1
487 END DO
488
489 IF (mod(counter, 2) == 0) THEN
490 determinant = z_one
491 ELSE
492 determinant = -z_one
493 END IF
494
495 ! compute product of diagonal elements
496 irow = 1
497 icol = 1
498 DO WHILE (irow <= nrow_local .AND. icol <= ncol_local)
499 IF (row_indices(irow) < col_indices(icol)) THEN
500 irow = irow + 1
501 ELSE IF (row_indices(irow) > col_indices(icol)) THEN
502 icol = icol + 1
503 ELSE ! diagonal element
504 determinant = determinant*a(irow, icol)
505 irow = irow + 1
506 icol = icol + 1
507 END IF
508 END DO
509 CALL matrix_a%matrix_struct%para_env%prod(determinant)
510#else
511 lda = SIZE(a, 1)
512 CALL zgetrf(nrow_global, nrow_global, a(1, 1), lda, ipivot, info)
513 counter = 0
514 determinant = z_one
515 DO irow = 1, nrow_global
516 IF (ipivot(irow) /= irow) counter = counter + 1
517 determinant = determinant*a(irow, irow)
518 END DO
519 IF (mod(counter, 2) == 1) determinant = -1.0_dp*determinant
520#endif
521
522 ! info is allowed to be zero
523 ! this does just signal a zero diagonal element
524 DEALLOCATE (ipivot)
525
526 CALL timestop(handle)
527 END SUBROUTINE cp_cfm_lu_decompose
528
529! **************************************************************************************************
530!> \brief Performs one of the matrix-matrix operations:
531!> matrix_c = alpha * op1( matrix_a ) * op2( matrix_b ) + beta*matrix_c.
532!> \param transa form of op1( matrix_a ):
533!> op1( matrix_a ) = matrix_a, when transa == 'N' ,
534!> op1( matrix_a ) = matrix_a^T, when transa == 'T' ,
535!> op1( matrix_a ) = matrix_a^H, when transa == 'C' ,
536!> \param transb form of op2( matrix_b )
537!> \param m number of rows of the matrix op1( matrix_a )
538!> \param n number of columns of the matrix op2( matrix_b )
539!> \param k number of columns of the matrix op1( matrix_a ) as well as
540!> number of rows of the matrix op2( matrix_b )
541!> \param alpha scale factor
542!> \param matrix_a matrix A
543!> \param matrix_b matrix B
544!> \param beta scale factor
545!> \param matrix_c matrix C
546!> \param a_first_col (optional) the first column of the matrix_a to multiply
547!> \param a_first_row (optional) the first row of the matrix_a to multiply
548!> \param b_first_col (optional) the first column of the matrix_b to multiply
549!> \param b_first_row (optional) the first row of the matrix_b to multiply
550!> \param c_first_col (optional) the first column of the matrix_c
551!> \param c_first_row (optional) the first row of the matrix_c
552!> \date 07.06.2001
553!> \author Matthias Krack
554!> \version 1.0
555! **************************************************************************************************
556 SUBROUTINE cp_cfm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, &
557 matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, &
558 c_first_row)
559 CHARACTER(len=1), INTENT(in) :: transa, transb
560 INTEGER, INTENT(in) :: m, n, k
561 COMPLEX(kind=dp), INTENT(in) :: alpha
562 TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, matrix_b
563 COMPLEX(kind=dp), INTENT(in) :: beta
564 TYPE(cp_cfm_type), INTENT(IN) :: matrix_c
565 INTEGER, INTENT(in), OPTIONAL :: a_first_col, a_first_row, b_first_col, &
566 b_first_row, c_first_col, c_first_row
567
568 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_gemm'
569
570 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b, c
571 INTEGER :: handle, i_a, i_b, i_c, j_a, j_b, j_c
572#if defined(__parallel)
573 INTEGER, DIMENSION(9) :: desca, descb, descc
574#else
575 INTEGER :: lda, ldb, ldc
576#endif
577
578 CALL timeset(routinen, handle)
579 a => matrix_a%local_data
580 b => matrix_b%local_data
581 c => matrix_c%local_data
582
583 i_a = 1
584 IF (PRESENT(a_first_row)) i_a = a_first_row
585
586 j_a = 1
587 IF (PRESENT(a_first_col)) j_a = a_first_col
588
589 i_b = 1
590 IF (PRESENT(b_first_row)) i_b = b_first_row
591
592 j_b = 1
593 IF (PRESENT(b_first_col)) j_b = b_first_col
594
595 i_c = 1
596 IF (PRESENT(c_first_row)) i_c = c_first_row
597
598 j_c = 1
599 IF (PRESENT(c_first_col)) j_c = c_first_col
600
601#if defined(__parallel)
602 desca(:) = matrix_a%matrix_struct%descriptor(:)
603 descb(:) = matrix_b%matrix_struct%descriptor(:)
604 descc(:) = matrix_c%matrix_struct%descriptor(:)
605
606 CALL pzgemm(transa, transb, m, n, k, alpha, a(1, 1), i_a, j_a, desca, &
607 b(1, 1), i_b, j_b, descb, beta, c(1, 1), i_c, j_c, descc)
608#else
609 lda = SIZE(a, 1)
610 ldb = SIZE(b, 1)
611 ldc = SIZE(c, 1)
612
613 ! consider zgemm3m
614 CALL zgemm(transa, transb, m, n, k, alpha, a(i_a, j_a), &
615 lda, b(i_b, j_b), ldb, beta, c(i_c, j_c), ldc)
616#endif
617 CALL timestop(handle)
618 END SUBROUTINE cp_cfm_gemm
619
620! **************************************************************************************************
621!> \brief Scales columns of the full matrix by corresponding factors.
622!> \param matrix_a matrix to scale
623!> \param scaling scale factors for every column. The actual number of scaled columns is
624!> limited by the number of scale factors given or by the actual number of columns
625!> whichever is smaller.
626!> \author Joost VandeVondele
627! **************************************************************************************************
628 SUBROUTINE cp_cfm_column_scale(matrix_a, scaling)
629 TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
630 COMPLEX(kind=dp), DIMENSION(:), INTENT(in) :: scaling
631
632 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_column_scale'
633
634 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
635 INTEGER :: handle, icol_local, ncol_local, &
636 nrow_local
637#if defined(__parallel)
638 INTEGER, DIMENSION(:), POINTER :: col_indices
639#endif
640
641 CALL timeset(routinen, handle)
642
643 a => matrix_a%local_data
644
645#if defined(__parallel)
646 CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, col_indices=col_indices)
647 ncol_local = min(ncol_local, SIZE(scaling))
648
649 DO icol_local = 1, ncol_local
650 CALL zscal(nrow_local, scaling(col_indices(icol_local)), a(1, icol_local), 1)
651 END DO
652#else
653 nrow_local = SIZE(a, 1)
654 ncol_local = min(SIZE(a, 2), SIZE(scaling))
655
656 DO icol_local = 1, ncol_local
657 CALL zscal(nrow_local, scaling(icol_local), a(1, icol_local), 1)
658 END DO
659#endif
660
661 CALL timestop(handle)
662 END SUBROUTINE cp_cfm_column_scale
663
664! **************************************************************************************************
665!> \brief Scales a complex matrix by a real number.
666!> matrix_a = alpha * matrix_b
667!> \param alpha scale factor
668!> \param matrix_a complex matrix to scale
669! **************************************************************************************************
670 SUBROUTINE cp_cfm_dscale(alpha, matrix_a)
671 REAL(kind=dp), INTENT(in) :: alpha
672 TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
673
674 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_dscale'
675
676 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
677 INTEGER :: handle
678
679 CALL timeset(routinen, handle)
680
681 NULLIFY (a)
682
683 a => matrix_a%local_data
684
685 CALL zdscal(SIZE(a), alpha, a(1, 1), 1)
686
687 CALL timestop(handle)
688 END SUBROUTINE cp_cfm_dscale
689
690! **************************************************************************************************
691!> \brief Scales a complex matrix by a complex number.
692!> matrix_a = alpha * matrix_b
693!> \param alpha scale factor
694!> \param matrix_a complex matrix to scale
695!> \note
696!> use cp_fm_set_all to zero (avoids problems with nan)
697! **************************************************************************************************
698 SUBROUTINE cp_cfm_zscale(alpha, matrix_a)
699 COMPLEX(kind=dp), INTENT(IN) :: alpha
700 TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
701
702 CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_zscale'
703
704 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
705 INTEGER :: handle, size_a
706
707 CALL timeset(routinen, handle)
708
709 NULLIFY (a)
710
711 a => matrix_a%local_data
712 size_a = SIZE(a, 1)*SIZE(a, 2)
713
714 CALL zscal(size_a, alpha, a(1, 1), 1)
715
716 CALL timestop(handle)
717 END SUBROUTINE cp_cfm_zscale
718
719! **************************************************************************************************
720!> \brief Solve the system of linear equations A*b=A_general using LU decomposition.
721!> Pay attention that both matrices are overwritten on exit and that
722!> the result is stored into the matrix 'general_a'.
723!> \param matrix_a matrix A (overwritten on exit)
724!> \param general_a (input) matrix A_general, (output) matrix B
725!> \param determinant (optional) determinant
726!> \author Florian Schiffmann
727! **************************************************************************************************
728 SUBROUTINE cp_cfm_solve(matrix_a, general_a, determinant)
729 TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, general_a
730 COMPLEX(kind=dp), OPTIONAL :: determinant
731
732 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_solve'
733
734 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, a_general
735 INTEGER :: counter, handle, info, irow, nrow_global
736 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
737
738#if defined(__parallel)
739 INTEGER :: icol, ncol_local, nrow_local
740 INTEGER, DIMENSION(9) :: desca, descb
741 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
742#else
743 INTEGER :: lda, ldb
744#endif
745
746 CALL timeset(routinen, handle)
747
748 a => matrix_a%local_data
749 a_general => general_a%local_data
750 nrow_global = matrix_a%matrix_struct%nrow_global
751 ALLOCATE (ipivot(nrow_global))
752
753#if defined(__parallel)
754 desca(:) = matrix_a%matrix_struct%descriptor(:)
755 descb(:) = general_a%matrix_struct%descriptor(:)
756 CALL pzgetrf(nrow_global, nrow_global, a(1, 1), 1, 1, desca, ipivot, info)
757 IF (PRESENT(determinant)) THEN
758 CALL cp_cfm_get_info(matrix_a, nrow_local=nrow_local, ncol_local=ncol_local, &
759 row_indices=row_indices, col_indices=col_indices)
760
761 counter = 0
762 DO irow = 1, nrow_local
763 IF (ipivot(irow) /= row_indices(irow)) counter = counter + 1
764 END DO
765
766 IF (mod(counter, 2) == 0) THEN
767 determinant = z_one
768 ELSE
769 determinant = -z_one
770 END IF
771
772 ! compute product of diagonal elements
773 irow = 1
774 icol = 1
775 DO WHILE (irow <= nrow_local .AND. icol <= ncol_local)
776 IF (row_indices(irow) < col_indices(icol)) THEN
777 irow = irow + 1
778 ELSE IF (row_indices(irow) > col_indices(icol)) THEN
779 icol = icol + 1
780 ELSE ! diagonal element
781 determinant = determinant*a(irow, icol)
782 irow = irow + 1
783 icol = icol + 1
784 END IF
785 END DO
786 CALL matrix_a%matrix_struct%para_env%prod(determinant)
787 END IF
788
789 CALL pzgetrs("N", nrow_global, nrow_global, a(1, 1), 1, 1, desca, &
790 ipivot, a_general(1, 1), 1, 1, descb, info)
791#else
792 lda = SIZE(a, 1)
793 ldb = SIZE(a_general, 1)
794 CALL zgetrf(nrow_global, nrow_global, a(1, 1), lda, ipivot, info)
795 IF (PRESENT(determinant)) THEN
796 counter = 0
797 determinant = z_one
798 DO irow = 1, nrow_global
799 IF (ipivot(irow) /= irow) counter = counter + 1
800 determinant = determinant*a(irow, irow)
801 END DO
802 IF (mod(counter, 2) == 1) determinant = -1.0_dp*determinant
803 END IF
804 CALL zgetrs("N", nrow_global, nrow_global, a(1, 1), lda, ipivot, a_general(1, 1), ldb, info)
805#endif
806
807 ! info is allowed to be zero
808 ! this does just signal a zero diagonal element
809 DEALLOCATE (ipivot)
810 CALL timestop(handle)
811
812 END SUBROUTINE cp_cfm_solve
813
814! **************************************************************************************************
815!> \brief Inverts a matrix using LU decomposition. The input matrix will be overwritten.
816!> \param matrix input a general square non-singular matrix, outputs its inverse
817!> \param info_out optional, if present outputs the info from (p)zgetri
818!> \author Lianheng Tong
819! **************************************************************************************************
820 SUBROUTINE cp_cfm_lu_invert(matrix, info_out)
821 TYPE(cp_cfm_type), INTENT(IN) :: matrix
822 INTEGER, INTENT(out), OPTIONAL :: info_out
823
824 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_lu_invert'
825
826 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
827 COMPLEX(kind=dp), DIMENSION(1) :: work1
828 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: mat
829 INTEGER :: handle, info, lwork, nrows_global
830 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
831
832#if defined(__parallel)
833 INTEGER :: liwork
834 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
835 INTEGER, DIMENSION(1) :: iwork1
836 INTEGER, DIMENSION(9) :: desca
837#else
838 INTEGER :: lda
839#endif
840
841 CALL timeset(routinen, handle)
842
843 mat => matrix%local_data
844 nrows_global = matrix%matrix_struct%nrow_global
845 cpassert(nrows_global == matrix%matrix_struct%ncol_global)
846 ALLOCATE (ipivot(nrows_global))
847
848 ! do LU decomposition
849#if defined(__parallel)
850 desca = matrix%matrix_struct%descriptor
851 CALL pzgetrf(nrows_global, nrows_global, &
852 mat(1, 1), 1, 1, desca, ipivot, info)
853#else
854 lda = SIZE(mat, 1)
855 CALL zgetrf(nrows_global, nrows_global, &
856 mat(1, 1), lda, ipivot, info)
857#endif
858 IF (info /= 0) THEN
859 CALL cp_abort(__location__, "LU decomposition has failed")
860 END IF
861
862 ! do inversion
863#if defined(__parallel)
864 CALL pzgetri(nrows_global, mat(1, 1), 1, 1, desca, &
865 ipivot, work1, -1, iwork1, -1, info)
866 lwork = int(work1(1))
867 liwork = int(iwork1(1))
868 ALLOCATE (work(lwork))
869 ALLOCATE (iwork(liwork))
870 CALL pzgetri(nrows_global, mat(1, 1), 1, 1, desca, &
871 ipivot, work, lwork, iwork, liwork, info)
872 DEALLOCATE (iwork)
873#else
874 CALL zgetri(nrows_global, mat(1, 1), lda, ipivot, work1, -1, info)
875 lwork = int(work1(1))
876 ALLOCATE (work(lwork))
877 CALL zgetri(nrows_global, mat(1, 1), lda, ipivot, work, lwork, info)
878#endif
879 DEALLOCATE (work)
880 DEALLOCATE (ipivot)
881
882 IF (PRESENT(info_out)) THEN
883 info_out = info
884 ELSE
885 IF (info /= 0) &
886 CALL cp_abort(__location__, "LU inversion has failed")
887 END IF
888
889 CALL timestop(handle)
890
891 END SUBROUTINE cp_cfm_lu_invert
892
893! **************************************************************************************************
894!> \brief Returns the trace of matrix_a^T matrix_b, i.e
895!> sum_{i,j}(matrix_a(i,j)*matrix_b(i,j)) .
896!> \param matrix_a a complex matrix
897!> \param matrix_b another complex matrix
898!> \param trace value of the trace operator
899!> \par History
900!> * 09.2017 created [Sergey Chulkov]
901!> \author Sergey Chulkov
902!> \note
903!> Based on the subroutine cp_fm_trace(). Note the transposition of matrix_a!
904! **************************************************************************************************
905 SUBROUTINE cp_cfm_trace(matrix_a, matrix_b, trace)
906 TYPE(cp_cfm_type), INTENT(IN) :: matrix_a, matrix_b
907 COMPLEX(kind=dp), INTENT(out) :: trace
908
909 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_trace'
910
911 INTEGER :: handle, mypcol, myprow, ncol_local, &
912 npcol, nprow, nrow_local
913 TYPE(cp_blacs_env_type), POINTER :: context
914 TYPE(mp_comm_type) :: group
915
916 CALL timeset(routinen, handle)
917
918 context => matrix_a%matrix_struct%context
919 myprow = context%mepos(1)
920 mypcol = context%mepos(2)
921 nprow = context%num_pe(1)
922 npcol = context%num_pe(2)
923
924 group = matrix_a%matrix_struct%para_env
925
926 nrow_local = min(matrix_a%matrix_struct%nrow_locals(myprow), matrix_b%matrix_struct%nrow_locals(myprow))
927 ncol_local = min(matrix_a%matrix_struct%ncol_locals(mypcol), matrix_b%matrix_struct%ncol_locals(mypcol))
928
929 ! compute an accurate dot-product
930 trace = accurate_dot_product(matrix_a%local_data(1:nrow_local, 1:ncol_local), &
931 matrix_b%local_data(1:nrow_local, 1:ncol_local))
932
933 CALL group%sum(trace)
934
935 CALL timestop(handle)
936
937 END SUBROUTINE cp_cfm_trace
938
939! **************************************************************************************************
940!> \brief Multiplies in place by a triangular matrix:
941!> matrix_b = alpha op(triangular_matrix) matrix_b
942!> or (if side='R')
943!> matrix_b = alpha matrix_b op(triangular_matrix)
944!> op(triangular_matrix) is:
945!> triangular_matrix (if transa="N" and invert_tr=.false.)
946!> triangular_matrix^T (if transa="T" and invert_tr=.false.)
947!> triangular_matrix^H (if transa="C" and invert_tr=.false.)
948!> triangular_matrix^(-1) (if transa="N" and invert_tr=.true.)
949!> triangular_matrix^(-T) (if transa="T" and invert_tr=.true.)
950!> triangular_matrix^(-H) (if transa="C" and invert_tr=.true.)
951!> \param triangular_matrix the triangular matrix that multiplies the other
952!> \param matrix_b the matrix that gets multiplied and stores the result
953!> \param side on which side of matrix_b stays op(triangular_matrix)
954!> (defaults to 'L')
955!> \param transa_tr ...
956!> \param invert_tr if the triangular matrix should be inverted
957!> (defaults to false)
958!> \param uplo_tr if triangular_matrix is stored in the upper ('U') or
959!> lower ('L') triangle (defaults to 'U')
960!> \param unit_diag_tr if the diagonal elements of triangular_matrix should
961!> be assumed to be 1 (defaults to false)
962!> \param n_rows the number of rows of the result (defaults to
963!> size(matrix_b,1))
964!> \param n_cols the number of columns of the result (defaults to
965!> size(matrix_b,2))
966!> \param alpha ...
967!> \par History
968!> 08.2002 created [fawzi]
969!> \author Fawzi Mohamed
970!> \note
971!> needs an mpi env
972! **************************************************************************************************
973 SUBROUTINE cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, &
974 transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, &
975 alpha)
976 TYPE(cp_cfm_type), INTENT(IN) :: triangular_matrix, matrix_b
977 CHARACTER, INTENT(in), OPTIONAL :: side, transa_tr
978 LOGICAL, INTENT(in), OPTIONAL :: invert_tr
979 CHARACTER, INTENT(in), OPTIONAL :: uplo_tr
980 LOGICAL, INTENT(in), OPTIONAL :: unit_diag_tr
981 INTEGER, INTENT(in), OPTIONAL :: n_rows, n_cols
982 COMPLEX(kind=dp), INTENT(in), OPTIONAL :: alpha
983
984 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_triangular_multiply'
985
986 CHARACTER :: side_char, transa, unit_diag, uplo
987 COMPLEX(kind=dp) :: al
988 INTEGER :: handle, m, n
989 LOGICAL :: invert
990
991 CALL timeset(routinen, handle)
992 side_char = 'L'
993 unit_diag = 'N'
994 uplo = 'U'
995 transa = 'N'
996 invert = .false.
997 al = z_one
998 CALL cp_cfm_get_info(matrix_b, nrow_global=m, ncol_global=n)
999 IF (PRESENT(side)) side_char = side
1000 IF (PRESENT(invert_tr)) invert = invert_tr
1001 IF (PRESENT(uplo_tr)) uplo = uplo_tr
1002 IF (PRESENT(unit_diag_tr)) THEN
1003 IF (unit_diag_tr) THEN
1004 unit_diag = 'U'
1005 ELSE
1006 unit_diag = 'N'
1007 END IF
1008 END IF
1009 IF (PRESENT(transa_tr)) transa = transa_tr
1010 IF (PRESENT(alpha)) al = alpha
1011 IF (PRESENT(n_rows)) m = n_rows
1012 IF (PRESENT(n_cols)) n = n_cols
1013
1014 IF (invert) THEN
1015
1016#if defined(__parallel)
1017 CALL pztrsm(side_char, uplo, transa, unit_diag, m, n, al, &
1018 triangular_matrix%local_data(1, 1), 1, 1, &
1019 triangular_matrix%matrix_struct%descriptor, &
1020 matrix_b%local_data(1, 1), 1, 1, &
1021 matrix_b%matrix_struct%descriptor(1))
1022#else
1023 CALL ztrsm(side_char, uplo, transa, unit_diag, m, n, al, &
1024 triangular_matrix%local_data(1, 1), &
1025 SIZE(triangular_matrix%local_data, 1), &
1026 matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1))
1027#endif
1028
1029 ELSE
1030
1031#if defined(__parallel)
1032 CALL pztrmm(side_char, uplo, transa, unit_diag, m, n, al, &
1033 triangular_matrix%local_data(1, 1), 1, 1, &
1034 triangular_matrix%matrix_struct%descriptor, &
1035 matrix_b%local_data(1, 1), 1, 1, &
1036 matrix_b%matrix_struct%descriptor(1))
1037#else
1038 CALL ztrmm(side_char, uplo, transa, unit_diag, m, n, al, &
1039 triangular_matrix%local_data(1, 1), &
1040 SIZE(triangular_matrix%local_data, 1), &
1041 matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1))
1042#endif
1043
1044 END IF
1045
1046 CALL timestop(handle)
1047
1048 END SUBROUTINE cp_cfm_triangular_multiply
1049
1050! **************************************************************************************************
1051!> \brief Inverts a triangular matrix.
1052!> \param matrix_a ...
1053!> \param uplo ...
1054!> \param info_out ...
1055!> \author MI
1056! **************************************************************************************************
1057 SUBROUTINE cp_cfm_triangular_invert(matrix_a, uplo, info_out)
1058 TYPE(cp_cfm_type), INTENT(IN) :: matrix_a
1059 CHARACTER, INTENT(in), OPTIONAL :: uplo
1060 INTEGER, INTENT(out), OPTIONAL :: info_out
1061
1062 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_triangular_invert'
1063
1064 CHARACTER :: unit_diag, my_uplo
1065 INTEGER :: handle, info, ncol_global
1066 COMPLEX(kind=dp), DIMENSION(:, :), &
1067 POINTER :: a
1068#if defined(__parallel)
1069 INTEGER, DIMENSION(9) :: desca
1070#endif
1071
1072 CALL timeset(routinen, handle)
1073
1074 unit_diag = 'N'
1075 my_uplo = 'U'
1076 IF (PRESENT(uplo)) my_uplo = uplo
1077
1078 ncol_global = matrix_a%matrix_struct%ncol_global
1079
1080 a => matrix_a%local_data
1081
1082#if defined(__parallel)
1083 desca(:) = matrix_a%matrix_struct%descriptor(:)
1084 CALL pztrtri(my_uplo, unit_diag, ncol_global, a(1, 1), 1, 1, desca, info)
1085#else
1086 CALL ztrtri(my_uplo, unit_diag, ncol_global, a(1, 1), ncol_global, info)
1087#endif
1088
1089 IF (PRESENT(info_out)) THEN
1090 info_out = info
1091 ELSE
1092 IF (info /= 0) &
1093 CALL cp_abort(__location__, &
1094 "triangular invert failed: matrix is not positive definite or ill-conditioned")
1095 END IF
1096
1097 CALL timestop(handle)
1098 END SUBROUTINE cp_cfm_triangular_invert
1099
1100! **************************************************************************************************
1101!> \brief Transposes a BLACS distributed complex matrix.
1102!> \param matrix input matrix
1103!> \param trans 'T' for transpose, 'C' for Hermitian conjugate
1104!> \param matrixt output matrix
1105!> \author Lianheng Tong
1106! **************************************************************************************************
1107 SUBROUTINE cp_cfm_transpose(matrix, trans, matrixt)
1108 TYPE(cp_cfm_type), INTENT(IN) :: matrix
1109 CHARACTER, INTENT(in) :: trans
1110 TYPE(cp_cfm_type), INTENT(IN) :: matrixt
1111
1112 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_transpose'
1113
1114 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa, cc
1115 INTEGER :: handle, ncol_global, nrow_global
1116#if defined(__parallel)
1117 INTEGER, DIMENSION(9) :: desca, descc
1118#elif !defined(__MKL)
1119 INTEGER :: ii, jj
1120#endif
1121
1122 CALL timeset(routinen, handle)
1123
1124 nrow_global = matrix%matrix_struct%nrow_global
1125 ncol_global = matrix%matrix_struct%ncol_global
1126
1127 cpassert(matrixt%matrix_struct%nrow_global == ncol_global)
1128 cpassert(matrixt%matrix_struct%ncol_global == nrow_global)
1129
1130 aa => matrix%local_data
1131 cc => matrixt%local_data
1132
1133#if defined(__parallel)
1134 desca = matrix%matrix_struct%descriptor
1135 descc = matrixt%matrix_struct%descriptor
1136 SELECT CASE (trans)
1137 CASE ('T')
1138 CALL pztranu(nrow_global, ncol_global, &
1139 z_one, aa(1, 1), 1, 1, desca, &
1140 z_zero, cc(1, 1), 1, 1, descc)
1141 CASE ('C')
1142 CALL pztranc(nrow_global, ncol_global, &
1143 z_one, aa(1, 1), 1, 1, desca, &
1144 z_zero, cc(1, 1), 1, 1, descc)
1145 CASE DEFAULT
1146 cpabort("trans only accepts 'T' or 'C'")
1147 END SELECT
1148#elif defined(__MKL)
1149 CALL mkl_zomatcopy('C', trans, nrow_global, ncol_global, 1.0_dp, aa(1, 1), nrow_global, cc(1, 1), ncol_global)
1150#else
1151 SELECT CASE (trans)
1152 CASE ('T')
1153 DO jj = 1, ncol_global
1154 DO ii = 1, nrow_global
1155 cc(ii, jj) = aa(jj, ii)
1156 END DO
1157 END DO
1158 CASE ('C')
1159 DO jj = 1, ncol_global
1160 DO ii = 1, nrow_global
1161 cc(ii, jj) = conjg(aa(jj, ii))
1162 END DO
1163 END DO
1164 CASE DEFAULT
1165 cpabort("trans only accepts 'T' or 'C'")
1166 END SELECT
1167#endif
1168
1169 CALL timestop(handle)
1170 END SUBROUTINE cp_cfm_transpose
1171
1172! **************************************************************************************************
1173!> \brief Norm of matrix using (p)zlange.
1174!> \param matrix input a general matrix
1175!> \param mode 'M' max abs element value,
1176!> '1' or 'O' one norm, i.e. maximum column sum,
1177!> 'I' infinity norm, i.e. maximum row sum,
1178!> 'F' or 'E' Frobenius norm, i.e. sqrt of sum of all squares of elements
1179!> \return the norm according to mode
1180!> \author Lianheng Tong
1181! **************************************************************************************************
1182 FUNCTION cp_cfm_norm(matrix, mode) RESULT(res)
1183 TYPE(cp_cfm_type), INTENT(IN) :: matrix
1184 CHARACTER, INTENT(IN) :: mode
1185 REAL(kind=dp) :: res
1186
1187 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_norm'
1188
1189 COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa
1190 INTEGER :: handle, lwork, ncols, ncols_local, &
1191 nrows, nrows_local
1192 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
1193
1194#if defined(__parallel)
1195 INTEGER, DIMENSION(9) :: desca
1196#else
1197 INTEGER :: lda
1198#endif
1199
1200 CALL timeset(routinen, handle)
1201
1202 CALL cp_cfm_get_info(matrix=matrix, &
1203 nrow_global=nrows, &
1204 ncol_global=ncols, &
1205 nrow_local=nrows_local, &
1206 ncol_local=ncols_local)
1207 aa => matrix%local_data
1208
1209 SELECT CASE (mode)
1210 CASE ('M', 'm')
1211 lwork = 1
1212 CASE ('1', 'O', 'o')
1213#if defined(__parallel)
1214 lwork = ncols_local
1215#else
1216 lwork = 1
1217#endif
1218 CASE ('I', 'i')
1219#if defined(__parallel)
1220 lwork = nrows_local
1221#else
1222 lwork = nrows
1223#endif
1224 CASE ('F', 'f', 'E', 'e')
1225 lwork = 1
1226 CASE DEFAULT
1227 cpabort("mode input is not valid")
1228 END SELECT
1229
1230 ALLOCATE (work(lwork))
1231
1232#if defined(__parallel)
1233 desca = matrix%matrix_struct%descriptor
1234 res = pzlange(mode, nrows, ncols, aa(1, 1), 1, 1, desca, work)
1235#else
1236 lda = SIZE(aa, 1)
1237 res = zlange(mode, nrows, ncols, aa(1, 1), lda, work)
1238#endif
1239
1240 DEALLOCATE (work)
1241 CALL timestop(handle)
1242 END FUNCTION cp_cfm_norm
1243
1244! **************************************************************************************************
1245!> \brief Applies a planar rotation defined by cs and sn to the i'th and j'th rows.
1246!> \param matrix ...
1247!> \param irow ...
1248!> \param jrow ...
1249!> \param cs cosine of the rotation angle
1250!> \param sn sinus of the rotation angle
1251!> \author Ole Schuett
1252! **************************************************************************************************
1253 SUBROUTINE cp_cfm_rot_rows(matrix, irow, jrow, cs, sn)
1254 TYPE(cp_cfm_type), INTENT(IN) :: matrix
1255 INTEGER, INTENT(IN) :: irow, jrow
1256 REAL(dp), INTENT(IN) :: cs, sn
1257
1258 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_rot_rows'
1259 INTEGER :: handle, ncol
1260 COMPLEX(KIND=dp) :: sn_cmplx
1261
1262#if defined(__parallel)
1263 INTEGER :: info, lwork
1264 INTEGER, DIMENSION(9) :: desc
1265 REAL(dp), DIMENSION(:), ALLOCATABLE :: work
1266#endif
1267 CALL timeset(routinen, handle)
1268 CALL cp_cfm_get_info(matrix, ncol_global=ncol)
1269 sn_cmplx = cmplx(sn, 0.0_dp, dp)
1270#if defined(__parallel)
1271 IF (1 /= matrix%matrix_struct%context%n_pid) THEN
1272 lwork = 2*ncol + 1
1273 ALLOCATE (work(lwork))
1274 desc(:) = matrix%matrix_struct%descriptor(:)
1275 info = 0
1276 CALL pzrot(ncol, &
1277 matrix%local_data(1, 1), irow, 1, desc, ncol, &
1278 matrix%local_data(1, 1), jrow, 1, desc, ncol, &
1279 cs, sn_cmplx, work, lwork, info)
1280 cpassert(info == 0)
1281 DEALLOCATE (work)
1282 ELSE
1283#endif
1284 CALL zrot(ncol, matrix%local_data(irow, 1), ncol, matrix%local_data(jrow, 1), ncol, cs, sn_cmplx)
1285#if defined(__parallel)
1286 END IF
1287#endif
1288 CALL timestop(handle)
1289 END SUBROUTINE cp_cfm_rot_rows
1290
1291! **************************************************************************************************
1292!> \brief Applies a planar rotation defined by cs and sn to the i'th and j'th columnns.
1293!> \param matrix ...
1294!> \param icol ...
1295!> \param jcol ...
1296!> \param cs cosine of the rotation angle
1297!> \param sn sinus of the rotation angle
1298!> \author Ole Schuett
1299! **************************************************************************************************
1300 SUBROUTINE cp_cfm_rot_cols(matrix, icol, jcol, cs, sn)
1301 TYPE(cp_cfm_type), INTENT(IN) :: matrix
1302 INTEGER, INTENT(IN) :: icol, jcol
1303 REAL(dp), INTENT(IN) :: cs, sn
1304
1305 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_rot_cols'
1306 INTEGER :: handle, nrow
1307 COMPLEX(KIND=dp) :: sn_cmplx
1308
1309#if defined(__parallel)
1310 INTEGER :: info, lwork
1311 INTEGER, DIMENSION(9) :: desc
1312 REAL(dp), DIMENSION(:), ALLOCATABLE :: work
1313#endif
1314 CALL timeset(routinen, handle)
1315 CALL cp_cfm_get_info(matrix, nrow_global=nrow)
1316 sn_cmplx = cmplx(sn, 0.0_dp, dp)
1317#if defined(__parallel)
1318 IF (1 /= matrix%matrix_struct%context%n_pid) THEN
1319 lwork = 2*nrow + 1
1320 ALLOCATE (work(lwork))
1321 desc(:) = matrix%matrix_struct%descriptor(:)
1322 info = 0
1323 CALL pzrot(nrow, &
1324 matrix%local_data(1, 1), 1, icol, desc, 1, &
1325 matrix%local_data(1, 1), 1, jcol, desc, 1, &
1326 cs, sn_cmplx, work, lwork, info)
1327 cpassert(info == 0)
1328 DEALLOCATE (work)
1329 ELSE
1330#endif
1331 CALL zrot(nrow, matrix%local_data(1, icol), 1, matrix%local_data(1, jcol), 1, cs, sn_cmplx)
1332#if defined(__parallel)
1333 END IF
1334#endif
1335 CALL timestop(handle)
1336 END SUBROUTINE cp_cfm_rot_cols
1337
1338! **************************************************************************************************
1339!> \brief ...
1340!> \param matrix ...
1341!> \param workspace ...
1342!> \param uplo triangular format; defaults to 'U'
1343!> \par History
1344!> 12.2024 Added optional workspace as input [Rocco Meli]
1345!> \author Jan Wilhelm
1346! **************************************************************************************************
1347 SUBROUTINE cp_cfm_uplo_to_full(matrix, workspace, uplo)
1348
1349 TYPE(cp_cfm_type), INTENT(IN) :: matrix
1350 TYPE(cp_cfm_type), INTENT(IN), OPTIONAL :: workspace
1351 CHARACTER, INTENT(IN), OPTIONAL :: uplo
1352
1353 CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_cfm_uplo_to_full'
1354
1355 CHARACTER :: myuplo
1356 INTEGER :: handle, i_global, iib, j_global, jjb, &
1357 ncol_local, nrow_local
1358 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1359 TYPE(cp_cfm_type) :: work
1360
1361 CALL timeset(routinen, handle)
1362
1363 IF (.NOT. PRESENT(workspace)) THEN
1364 CALL cp_cfm_create(work, matrix%matrix_struct)
1365 ELSE
1366 work = workspace
1367 END IF
1368
1369 myuplo = 'U'
1370 IF (PRESENT(uplo)) myuplo = uplo
1371
1372 ! get info of fm_mat_Q
1373 CALL cp_cfm_get_info(matrix=matrix, &
1374 nrow_local=nrow_local, &
1375 ncol_local=ncol_local, &
1376 row_indices=row_indices, &
1377 col_indices=col_indices)
1378
1379 DO jjb = 1, ncol_local
1380 j_global = col_indices(jjb)
1381 DO iib = 1, nrow_local
1382 i_global = row_indices(iib)
1383 IF (merge(j_global < i_global, j_global > i_global, (myuplo == "U") .OR. (myuplo == "u"))) THEN
1384 matrix%local_data(iib, jjb) = z_zero
1385 ELSE IF (j_global == i_global) THEN
1386 matrix%local_data(iib, jjb) = matrix%local_data(iib, jjb)/(2.0_dp, 0.0_dp)
1387 END IF
1388 END DO
1389 END DO
1390
1391 CALL cp_cfm_transpose(matrix, 'C', work)
1392
1393 CALL cp_cfm_scale_and_add(z_one, matrix, z_one, work)
1394
1395 IF (.NOT. PRESENT(workspace)) THEN
1396 CALL cp_cfm_release(work)
1397 END IF
1398
1399 CALL timestop(handle)
1400
1401 END SUBROUTINE cp_cfm_uplo_to_full
1402
1403END 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:56
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:821
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:1183
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:559
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:729
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:1108
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:1254
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:349
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:976
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:1348
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:629
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:1301
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:1058
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:454
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:906
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:182
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:125
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:630
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:69
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
message_passing::mp_comm_type
Definition message_passing.F:155