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