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