Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Basic linear algebra operations for full matrices.
10 : !> \par History
11 : !> 08.2002 split out of qs_blacs [fawzi]
12 : !> \author Fawzi Mohamed
13 : ! **************************************************************************************************
14 : MODULE cp_fm_basic_linalg
15 : USE cp_blacs_env, ONLY: cp_blacs_env_type
16 : USE cp_fm_struct, ONLY: cp_fm_struct_equivalent
17 : USE cp_fm_types, ONLY: &
18 : cp_fm_create, cp_fm_get_diag, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_p_type, &
19 : cp_fm_release, cp_fm_set_all, cp_fm_set_element, cp_fm_set_submatrix, cp_fm_to_fm, &
20 : cp_fm_type
21 : USE cp_log_handling, ONLY: cp_logger_get_default_unit_nr, &
22 : cp_to_string
23 : USE kahan_sum, ONLY: accurate_dot_product, &
24 : accurate_sum
25 : USE kinds, ONLY: dp, &
26 : int_8
27 : USE machine, ONLY: m_memory
28 : USE mathlib, ONLY: get_pseudo_inverse_svd, &
29 : invert_matrix
30 : USE message_passing, ONLY: mp_comm_type
31 : #if defined (__HAS_IEEE_EXCEPTIONS)
32 : USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
33 : ieee_set_halting_mode, &
34 : IEEE_ALL
35 : #endif
36 : #include "../base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 : PRIVATE
40 :
41 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_fm_basic_linalg'
43 :
44 : PUBLIC :: cp_fm_scale, & ! scale a matrix
45 : cp_fm_scale_and_add, & ! scale and add two matrices
46 : cp_fm_geadd, & ! general addition
47 : cp_fm_add_columns, & ! daxpy type column add
48 : cp_fm_column_scale, & ! scale columns of a matrix
49 : cp_fm_row_scale, & ! scale rows of a matrix
50 : cp_fm_trace, & ! trace of the transpose(A)*B
51 : cp_fm_contracted_trace, & ! sum_{i,...,k} Tr [A(i,...,k)^T * B(i,...,k)]
52 : cp_fm_norm, & ! different norms of A
53 : cp_fm_schur_product, & ! schur product
54 : cp_fm_transpose, & ! transpose a matrix
55 : cp_fm_uplo_to_full, & ! symmetrise a triangular matrix
56 : cp_fm_syrk, & ! rank k update
57 : cp_fm_triangular_multiply, & ! triangular matrix multiply / solve
58 : cp_fm_symm, & ! multiply a symmetric with a non-symmetric matrix
59 : cp_fm_gemm, & ! multiply two matrices
60 : cp_complex_fm_gemm, & ! multiply two complex matrices, represented by non_complex fm matrices
61 : cp_fm_invert, & ! computes the inverse and determinant
62 : cp_fm_frobenius_norm, & ! frobenius norm
63 : cp_fm_triangular_invert, & ! compute the reciprocal of a triangular matrix
64 : cp_fm_qr_factorization, & ! compute the QR factorization of a rectangular matrix
65 : cp_fm_solve, & ! solves the equation A*B=C A and C are input
66 : cp_fm_pdgeqpf, & ! compute a QR factorization with column pivoting of a M-by-N distributed matrix
67 : cp_fm_pdorgqr, & ! generates an M-by-N as first N columns of a product of K elementary reflectors
68 : cp_fm_potrf, & ! Cholesky decomposition
69 : cp_fm_potri, & ! Invert triangular matrix
70 : cp_fm_rot_rows, & ! rotates two rows
71 : cp_fm_rot_cols, & ! rotates two columns
72 : cp_fm_Gram_Schmidt_orthonorm, & ! Gram-Schmidt orthonormalization of columns of a full matrix, &
73 : cp_fm_det, & ! determinant of a real matrix with correct sign
74 : cp_fm_matvec, & ! matrix-vector multiplication (vector replicated)
75 : cp_fm_add_on_diag ! add a scalar to the diagonal
76 :
77 : REAL(KIND=dp), EXTERNAL :: dlange, pdlange, pdlatra
78 :
79 : INTERFACE cp_fm_trace
80 : MODULE PROCEDURE cp_fm_trace_a0b0t0
81 : MODULE PROCEDURE cp_fm_trace_a1b0t1_a
82 : MODULE PROCEDURE cp_fm_trace_a1b0t1_p
83 : MODULE PROCEDURE cp_fm_trace_a1b1t1_aa
84 : MODULE PROCEDURE cp_fm_trace_a1b1t1_ap
85 : MODULE PROCEDURE cp_fm_trace_a1b1t1_pa
86 : MODULE PROCEDURE cp_fm_trace_a1b1t1_pp
87 : END INTERFACE cp_fm_trace
88 :
89 : INTERFACE cp_fm_contracted_trace
90 : MODULE PROCEDURE cp_fm_contracted_trace_a2b2t2_aa
91 : MODULE PROCEDURE cp_fm_contracted_trace_a2b2t2_ap
92 : MODULE PROCEDURE cp_fm_contracted_trace_a2b2t2_pa
93 : MODULE PROCEDURE cp_fm_contracted_trace_a2b2t2_pp
94 : END INTERFACE cp_fm_contracted_trace
95 : CONTAINS
96 :
97 : ! **************************************************************************************************
98 : !> \brief Computes the determinant (with a correct sign even in parallel environment!) of a real square matrix
99 : !> \author A. Sinyavskiy (andrey.sinyavskiy@chem.uzh.ch)
100 : ! **************************************************************************************************
101 0 : SUBROUTINE cp_fm_det(matrix_a, det_a)
102 :
103 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a
104 : REAL(KIND=dp), INTENT(OUT) :: det_a
105 : REAL(KIND=dp) :: determinant
106 : TYPE(cp_fm_type) :: matrix_lu
107 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
108 : INTEGER :: n, i, info, P
109 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
110 : REAL(KIND=dp), DIMENSION(:), POINTER :: diag
111 :
112 : #if defined(__parallel)
113 : INTEGER :: myprow, nprow, npcol, nrow_local, nrow_block, irow_local
114 : INTEGER, DIMENSION(9) :: desca
115 : #endif
116 :
117 : CALL cp_fm_create(matrix=matrix_lu, &
118 : matrix_struct=matrix_a%matrix_struct, &
119 0 : name="A_lu"//TRIM(ADJUSTL(cp_to_string(1)))//"MATRIX")
120 0 : CALL cp_fm_to_fm(matrix_a, matrix_lu)
121 :
122 0 : a => matrix_lu%local_data
123 0 : n = matrix_lu%matrix_struct%nrow_global
124 0 : ALLOCATE (ipivot(n))
125 0 : ipivot(:) = 0
126 0 : P = 0
127 0 : ALLOCATE (diag(n))
128 0 : diag(:) = 0.0_dp
129 : #if defined(__parallel)
130 : ! Use LU decomposition
131 0 : desca(:) = matrix_lu%matrix_struct%descriptor(:)
132 0 : CALL pdgetrf(n, n, a, 1, 1, desca, ipivot, info)
133 0 : CALL cp_fm_get_diag(matrix_lu, diag)
134 0 : determinant = PRODUCT(diag)
135 0 : myprow = matrix_lu%matrix_struct%context%mepos(1)
136 0 : nprow = matrix_lu%matrix_struct%context%num_pe(1)
137 0 : npcol = matrix_lu%matrix_struct%context%num_pe(2)
138 0 : nrow_local = matrix_lu%matrix_struct%nrow_locals(myprow)
139 0 : nrow_block = matrix_lu%matrix_struct%nrow_block
140 0 : DO irow_local = 1, nrow_local
141 0 : i = matrix_lu%matrix_struct%row_indices(irow_local)
142 0 : IF (ipivot(irow_local) /= i) P = P + 1
143 : END DO
144 0 : CALL matrix_lu%matrix_struct%para_env%sum(P)
145 : ! very important fix
146 0 : P = P/npcol
147 : #else
148 : CALL dgetrf(n, n, a, n, ipivot, info)
149 : CALL cp_fm_get_diag(matrix_lu, diag)
150 : determinant = PRODUCT(diag)
151 : DO i = 1, n
152 : IF (ipivot(i) /= i) P = P + 1
153 : END DO
154 : #endif
155 0 : DEALLOCATE (ipivot)
156 0 : DEALLOCATE (diag)
157 0 : CALL cp_fm_release(matrix_lu)
158 0 : det_a = determinant*(-2*MOD(P, 2) + 1.0_dp)
159 0 : END SUBROUTINE cp_fm_det
160 :
161 : ! **************************************************************************************************
162 : !> \brief calc A <- alpha*A + beta*B
163 : !> optimized for alpha == 1.0 (just add beta*B) and beta == 0.0 (just
164 : !> scale A)
165 : !> \param alpha ...
166 : !> \param matrix_a ...
167 : !> \param beta ...
168 : !> \param matrix_b ...
169 : ! **************************************************************************************************
170 1397586 : SUBROUTINE cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
171 :
172 : REAL(KIND=dp), INTENT(IN) :: alpha
173 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a
174 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: beta
175 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: matrix_b
176 :
177 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_scale_and_add'
178 :
179 : INTEGER :: handle, size_a, size_b
180 : REAL(KIND=dp) :: my_beta
181 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
182 :
183 1397586 : CALL timeset(routineN, handle)
184 :
185 1397586 : my_beta = 0.0_dp
186 1397586 : IF (PRESENT(matrix_b)) my_beta = 1.0_dp
187 1397586 : IF (PRESENT(beta)) my_beta = beta
188 : NULLIFY (a, b)
189 :
190 1397586 : IF (PRESENT(beta)) THEN
191 1397584 : CPASSERT(PRESENT(matrix_b))
192 1397584 : IF (ASSOCIATED(matrix_a%local_data, matrix_b%local_data)) THEN
193 0 : CPWARN("Bad use of routine. Call cp_fm_scale instead")
194 0 : CALL cp_fm_scale(alpha + beta, matrix_a)
195 0 : CALL timestop(handle)
196 0 : RETURN
197 : END IF
198 : END IF
199 :
200 1397586 : a => matrix_a%local_data
201 :
202 1397586 : size_a = SIZE(a, 1)*SIZE(a, 2)
203 :
204 1397586 : IF (alpha /= 1.0_dp) THEN
205 85308 : CALL dscal(size_a, alpha, a, 1)
206 : END IF
207 1397586 : IF (my_beta /= 0.0_dp) THEN
208 1389650 : IF (matrix_a%matrix_struct%context /= matrix_b%matrix_struct%context) &
209 0 : CPABORT("Matrices must be in the same blacs context")
210 :
211 1389650 : IF (cp_fm_struct_equivalent(matrix_a%matrix_struct, &
212 : matrix_b%matrix_struct)) THEN
213 :
214 1389650 : b => matrix_b%local_data
215 1389650 : size_b = SIZE(b, 1)*SIZE(b, 2)
216 1389650 : IF (size_a /= size_b) &
217 0 : CPABORT("Matrices must have same local sizes")
218 :
219 1389650 : CALL daxpy(size_a, my_beta, b, 1, a, 1)
220 :
221 : ELSE
222 : CALL cp_abort(__LOCATION__, &
223 : "cp_fm_scale_and_add is not yet implemented for cases "// &
224 0 : "where two input matrix structures are not equivalent")
225 : END IF
226 :
227 : END IF
228 :
229 1397586 : CALL timestop(handle)
230 :
231 : END SUBROUTINE cp_fm_scale_and_add
232 :
233 : ! **************************************************************************************************
234 : !> \brief interface to BLACS geadd:
235 : !> matrix_b = beta*matrix_b + alpha*opt(matrix_a)
236 : !> where opt(matrix_a) can be either:
237 : !> 'N': matrix_a
238 : !> 'T': matrix_a^T
239 : !> 'C': matrix_a^H (Hermitian conjugate)
240 : !> note that this is a level three routine, use cp_fm_scale_and_add if that
241 : !> is sufficient for your needs
242 : !> \param alpha : complex scalar
243 : !> \param trans : 'N' normal, 'T' transposed
244 : !> \param matrix_a : input matrix_a
245 : !> \param beta : complex scalar
246 : !> \param matrix_b : input matrix_b, upon out put the updated matrix_b
247 : !> \author Lianheng Tong
248 : ! **************************************************************************************************
249 650 : SUBROUTINE cp_fm_geadd(alpha, trans, matrix_a, beta, matrix_b)
250 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
251 : CHARACTER, INTENT(IN) :: trans
252 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a, matrix_b
253 :
254 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_geadd'
255 :
256 : INTEGER :: nrow_global, ncol_global, handle
257 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: aa, bb
258 : #if defined(__parallel)
259 : INTEGER, DIMENSION(9) :: desca, descb
260 : #elif !defined(__MKL)
261 : INTEGER :: ii, jj
262 : #endif
263 :
264 650 : CALL timeset(routineN, handle)
265 :
266 650 : nrow_global = matrix_a%matrix_struct%nrow_global
267 650 : ncol_global = matrix_a%matrix_struct%ncol_global
268 650 : CPASSERT(nrow_global == matrix_b%matrix_struct%nrow_global)
269 650 : CPASSERT(ncol_global == matrix_b%matrix_struct%ncol_global)
270 :
271 650 : aa => matrix_a%local_data
272 650 : bb => matrix_b%local_data
273 :
274 : #if defined(__parallel)
275 6500 : desca = matrix_a%matrix_struct%descriptor
276 6500 : descb = matrix_b%matrix_struct%descriptor
277 : CALL pdgeadd(trans, &
278 : nrow_global, &
279 : ncol_global, &
280 : alpha, &
281 : aa, &
282 : 1, 1, &
283 : desca, &
284 : beta, &
285 : bb, &
286 : 1, 1, &
287 650 : descb)
288 : #elif defined(__MKL)
289 : CALL mkl_domatadd('C', trans, 'N', nrow_global, ncol_global, &
290 : alpha, aa, nrow_global, beta, bb, nrow_global, bb, nrow_global)
291 : #else
292 : ! dgeadd is not a standard BLAS function, although it is implemented
293 : ! in some libraries like OpenBLAS, so not going to use it here
294 : SELECT CASE (trans)
295 : CASE ('T')
296 : DO jj = 1, ncol_global
297 : DO ii = 1, nrow_global
298 : bb(ii, jj) = beta*bb(ii, jj) + alpha*aa(jj, ii)
299 : END DO
300 : END DO
301 : CASE DEFAULT
302 : DO jj = 1, ncol_global
303 : DO ii = 1, nrow_global
304 : bb(ii, jj) = beta*bb(ii, jj) + alpha*aa(ii, jj)
305 : END DO
306 : END DO
307 : END SELECT
308 : #endif
309 :
310 650 : CALL timestop(handle)
311 :
312 650 : END SUBROUTINE cp_fm_geadd
313 :
314 : ! **************************************************************************************************
315 : !> \brief Add (and scale) a subset of columns of a fm to a fm
316 : !> b = alpha*a + b
317 : !> \param msource ...
318 : !> \param mtarget ...
319 : !> \param ncol ...
320 : !> \param alpha ...
321 : !> \param source_start ...
322 : !> \param target_start ...
323 : ! **************************************************************************************************
324 0 : SUBROUTINE cp_fm_add_columns(msource, mtarget, ncol, alpha, source_start, target_start)
325 :
326 : TYPE(cp_fm_type), INTENT(IN) :: msource, mtarget
327 : INTEGER, INTENT(IN) :: ncol
328 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: alpha
329 : INTEGER, INTENT(IN), OPTIONAL :: source_start, target_start
330 :
331 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_add_columns'
332 :
333 : INTEGER :: handle, n, ss, ts, i
334 : REAL(KIND=dp) :: fscale
335 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
336 : #if defined(__parallel)
337 : INTEGER, DIMENSION(9) :: desca, descb
338 : #endif
339 :
340 0 : CALL timeset(routineN, handle)
341 :
342 0 : ss = 1
343 0 : ts = 1
344 :
345 0 : IF (PRESENT(source_start)) ss = source_start
346 0 : IF (PRESENT(target_start)) ts = target_start
347 :
348 0 : fscale = 1.0_dp
349 0 : IF (PRESENT(alpha)) fscale = alpha
350 :
351 0 : n = msource%matrix_struct%nrow_global
352 :
353 0 : a => msource%local_data
354 0 : b => mtarget%local_data
355 :
356 : #if defined(__parallel)
357 : MARK_USED(i)
358 0 : desca(:) = msource%matrix_struct%descriptor(:)
359 0 : descb(:) = mtarget%matrix_struct%descriptor(:)
360 0 : CALL pdgeadd("N", n, ncol, fscale, a, 1, ss, desca, 1.0_dp, b, 1, ts, descb)
361 : #else
362 : DO i = 0, ncol - 1
363 : b(1:n, ts + i) = b(1:n, ts + i) + fscale*a(1:n, ss + i)
364 : END DO
365 : #endif
366 :
367 0 : CALL timestop(handle)
368 :
369 0 : END SUBROUTINE cp_fm_add_columns
370 :
371 : ! **************************************************************************************************
372 : !> \brief Computes the LU-decomposition of the matrix, and the determinant of the matrix
373 : !> IMPORTANT : the sign of the determinant is not defined correctly yet ....
374 : !> \param matrix_a ...
375 : !> \param almost_determinant ...
376 : !> \param correct_sign ...
377 : !> \par History
378 : !> added correct_sign 02.07 (fschiff)
379 : !> \author Joost VandeVondele
380 : !> \note
381 : !> - matrix_a is overwritten
382 : !> - the sign of the determinant might be wrong
383 : !> - SERIOUS WARNING (KNOWN BUG) : the sign of the determinant depends on ipivot
384 : !> - one should be able to find out if ipivot is an even or an odd permutation...
385 : !> if you need the correct sign, just add correct_sign==.TRUE. (fschiff)
386 : !> - Use cp_fm_get_diag instead of n times cp_fm_get_element (A. Bussy)
387 : ! **************************************************************************************************
388 0 : SUBROUTINE cp_fm_lu_decompose(matrix_a, almost_determinant, correct_sign)
389 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a
390 : REAL(KIND=dp), INTENT(OUT) :: almost_determinant
391 : LOGICAL, INTENT(IN), OPTIONAL :: correct_sign
392 :
393 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_lu_decompose'
394 :
395 : INTEGER :: handle, i, info, n
396 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
397 : REAL(KIND=dp) :: determinant
398 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
399 : #if defined(__parallel)
400 : INTEGER, DIMENSION(9) :: desca
401 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: diag
402 : #else
403 : INTEGER :: lda
404 : #endif
405 :
406 0 : CALL timeset(routineN, handle)
407 :
408 0 : a => matrix_a%local_data
409 0 : n = matrix_a%matrix_struct%nrow_global
410 0 : ALLOCATE (ipivot(n + matrix_a%matrix_struct%nrow_block))
411 :
412 : #if defined(__parallel)
413 : MARK_USED(correct_sign)
414 0 : desca(:) = matrix_a%matrix_struct%descriptor(:)
415 0 : CALL pdgetrf(n, n, a, 1, 1, desca, ipivot, info)
416 :
417 0 : ALLOCATE (diag(n))
418 0 : diag(:) = 0.0_dp
419 0 : CALL cp_fm_get_diag(matrix_a, diag)
420 0 : determinant = 1.0_dp
421 0 : DO i = 1, n
422 0 : determinant = determinant*diag(i)
423 : END DO
424 0 : DEALLOCATE (diag)
425 : #else
426 : lda = SIZE(a, 1)
427 : CALL dgetrf(n, n, a, lda, ipivot, info)
428 : determinant = 1.0_dp
429 : IF (correct_sign) THEN
430 : DO i = 1, n
431 : IF (ipivot(i) /= i) THEN
432 : determinant = -determinant*a(i, i)
433 : ELSE
434 : determinant = determinant*a(i, i)
435 : END IF
436 : END DO
437 : ELSE
438 : DO i = 1, n
439 : determinant = determinant*a(i, i)
440 : END DO
441 : END IF
442 : #endif
443 : ! info is allowed to be zero
444 : ! this does just signal a zero diagonal element
445 0 : DEALLOCATE (ipivot)
446 0 : almost_determinant = determinant ! notice that the sign is random
447 0 : CALL timestop(handle)
448 0 : END SUBROUTINE cp_fm_lu_decompose
449 :
450 : ! **************************************************************************************************
451 : !> \brief computes matrix_c = beta * matrix_c + alpha * ( matrix_a ** transa ) * ( matrix_b ** transb )
452 : !> \param transa : 'N' -> normal 'T' -> transpose
453 : !> alpha,beta :: can be 0.0_dp and 1.0_dp
454 : !> \param transb same operation choices as transa
455 : !> \param m ...
456 : !> \param n ...
457 : !> \param k ...
458 : !> \param alpha ...
459 : !> \param matrix_a : m x k matrix ( ! for transa = 'N')
460 : !> \param matrix_b : k x n matrix ( ! for transb = 'N')
461 : !> \param beta ...
462 : !> \param matrix_c : m x n matrix
463 : !> \param a_first_col ...
464 : !> \param a_first_row ...
465 : !> \param b_first_col : the k x n matrix starts at col b_first_col of matrix_b (avoid usage)
466 : !> \param b_first_row ...
467 : !> \param c_first_col ...
468 : !> \param c_first_row ...
469 : !> \author Matthias Krack
470 : !> \note
471 : !> matrix_c should have no overlap with matrix_a, matrix_b
472 : ! **************************************************************************************************
473 19996 : SUBROUTINE cp_fm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, &
474 : matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, &
475 : c_first_col, c_first_row)
476 :
477 : CHARACTER(LEN=1), INTENT(IN) :: transa, transb
478 : INTEGER, INTENT(IN) :: m, n, k
479 : REAL(KIND=dp), INTENT(IN) :: alpha
480 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a, matrix_b
481 : REAL(KIND=dp), INTENT(IN) :: beta
482 : TYPE(cp_fm_type), INTENT(IN) :: matrix_c
483 : INTEGER, INTENT(IN), OPTIONAL :: a_first_col, a_first_row, &
484 : b_first_col, b_first_row, &
485 : c_first_col, c_first_row
486 :
487 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_gemm'
488 :
489 : INTEGER :: handle, i_a, i_b, i_c, j_a, &
490 : j_b, j_c
491 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, c
492 : #if defined(__parallel)
493 : INTEGER, DIMENSION(9) :: desca, descb, descc
494 : #else
495 : INTEGER :: lda, ldb, ldc
496 : #endif
497 :
498 19996 : CALL timeset(routineN, handle)
499 :
500 : !sample peak memory
501 19996 : CALL m_memory()
502 :
503 19996 : a => matrix_a%local_data
504 19996 : b => matrix_b%local_data
505 19996 : c => matrix_c%local_data
506 :
507 19996 : i_a = 1
508 19996 : IF (PRESENT(a_first_row)) i_a = a_first_row
509 :
510 19996 : j_a = 1
511 19996 : IF (PRESENT(a_first_col)) j_a = a_first_col
512 :
513 19996 : i_b = 1
514 19996 : IF (PRESENT(b_first_row)) i_b = b_first_row
515 :
516 19996 : j_b = 1
517 19996 : IF (PRESENT(b_first_col)) j_b = b_first_col
518 :
519 19996 : i_c = 1
520 19996 : IF (PRESENT(c_first_row)) i_c = c_first_row
521 :
522 19996 : j_c = 1
523 19996 : IF (PRESENT(c_first_col)) j_c = c_first_col
524 :
525 : #if defined(__parallel)
526 :
527 199960 : desca(:) = matrix_a%matrix_struct%descriptor(:)
528 199960 : descb(:) = matrix_b%matrix_struct%descriptor(:)
529 199960 : descc(:) = matrix_c%matrix_struct%descriptor(:)
530 :
531 : CALL pdgemm(transa, transb, m, n, k, alpha, a, i_a, j_a, desca, b, i_b, j_b, &
532 19996 : descb, beta, c, i_c, j_c, descc)
533 :
534 : #else
535 :
536 : lda = SIZE(a, 1)
537 : ldb = SIZE(b, 1)
538 : ldc = SIZE(c, 1)
539 :
540 : CALL dgemm(transa, transb, m, n, k, alpha, a(i_a, j_a), lda, b(i_b, j_b), ldb, beta, c(i_c, j_c), ldc)
541 :
542 : #endif
543 19996 : CALL timestop(handle)
544 :
545 19996 : END SUBROUTINE cp_fm_gemm
546 :
547 : ! **************************************************************************************************
548 : !> \brief computes matrix_c = beta * matrix_c + alpha * matrix_a * matrix_b
549 : !> computes matrix_c = beta * matrix_c + alpha * matrix_b * matrix_a
550 : !> where matrix_a is symmetric
551 : !> \param side : 'L' -> matrix_a is on the left 'R' -> matrix_a is on the right
552 : !> alpha,beta :: can be 0.0_dp and 1.0_dp
553 : !> \param uplo triangular format
554 : !> \param m ...
555 : !> \param n ...
556 : !> \param alpha ...
557 : !> \param matrix_a : m x m matrix
558 : !> \param matrix_b : m x n matrix
559 : !> \param beta ...
560 : !> \param matrix_c : m x n matrix
561 : !> \author Matthias Krack
562 : !> \note
563 : !> matrix_c should have no overlap with matrix_a, matrix_b
564 : !> all matrices in QS are triangular according to uplo
565 : !> matrix_a is always an m x m matrix
566 : !> typically slower than cp_fm_gemm (especially in parallel easily 50 percent)
567 : ! **************************************************************************************************
568 217224 : SUBROUTINE cp_fm_symm(side, uplo, m, n, alpha, matrix_a, matrix_b, beta, matrix_c)
569 :
570 : CHARACTER(LEN=1), INTENT(IN) :: side, uplo
571 : INTEGER, INTENT(IN) :: m, n
572 : REAL(KIND=dp), INTENT(IN) :: alpha
573 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a, matrix_b
574 : REAL(KIND=dp), INTENT(IN) :: beta
575 : TYPE(cp_fm_type), INTENT(IN) :: matrix_c
576 :
577 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_symm'
578 :
579 : INTEGER :: handle
580 217224 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, c
581 : #if defined(__parallel)
582 : INTEGER, DIMENSION(9) :: desca, descb, descc
583 : #else
584 : INTEGER :: lda, ldb, ldc
585 : #endif
586 :
587 217224 : CALL timeset(routineN, handle)
588 :
589 217224 : a => matrix_a%local_data
590 217224 : b => matrix_b%local_data
591 217224 : c => matrix_c%local_data
592 :
593 : #if defined(__parallel)
594 :
595 2172240 : desca(:) = matrix_a%matrix_struct%descriptor(:)
596 2172240 : descb(:) = matrix_b%matrix_struct%descriptor(:)
597 2172240 : descc(:) = matrix_c%matrix_struct%descriptor(:)
598 :
599 217224 : CALL pdsymm(side, uplo, m, n, alpha, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb, beta, c(1, 1), 1, 1, descc)
600 :
601 : #else
602 :
603 : lda = matrix_a%matrix_struct%local_leading_dimension
604 : ldb = matrix_b%matrix_struct%local_leading_dimension
605 : ldc = matrix_c%matrix_struct%local_leading_dimension
606 :
607 : CALL dsymm(side, uplo, m, n, alpha, a(1, 1), lda, b(1, 1), ldb, beta, c(1, 1), ldc)
608 :
609 : #endif
610 217224 : CALL timestop(handle)
611 :
612 217224 : END SUBROUTINE cp_fm_symm
613 :
614 : ! **************************************************************************************************
615 : !> \brief computes the Frobenius norm of matrix_a
616 : !> \brief computes the Frobenius norm of matrix_a
617 : !> \param matrix_a : m x n matrix
618 : !> \return ...
619 : !> \author VW
620 : ! **************************************************************************************************
621 7926 : FUNCTION cp_fm_frobenius_norm(matrix_a) RESULT(norm)
622 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a
623 : REAL(KIND=dp) :: norm
624 :
625 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_frobenius_norm'
626 :
627 : INTEGER :: handle, size_a
628 7926 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
629 : REAL(KIND=dp), EXTERNAL :: DDOT
630 : #if defined(__parallel)
631 : TYPE(mp_comm_type) :: group
632 : #endif
633 :
634 7926 : CALL timeset(routineN, handle)
635 :
636 : norm = 0.0_dp
637 7926 : a => matrix_a%local_data
638 7926 : size_a = SIZE(a, 1)*SIZE(a, 2)
639 7926 : norm = DDOT(size_a, a(1, 1), 1, a(1, 1), 1)
640 : #if defined(__parallel)
641 7926 : group = matrix_a%matrix_struct%para_env
642 7926 : CALL group%sum(norm)
643 : #endif
644 7926 : norm = SQRT(norm)
645 :
646 7926 : CALL timestop(handle)
647 :
648 7926 : END FUNCTION cp_fm_frobenius_norm
649 :
650 : ! **************************************************************************************************
651 : !> \brief performs a rank-k update of a symmetric matrix_c
652 : !> matrix_c = beta * matrix_c + alpha * matrix_a * transpose ( matrix_a )
653 : !> \param uplo : 'U' ('L')
654 : !> \param trans : 'N' ('T')
655 : !> \param k : number of cols to use in matrix_a
656 : !> ia,ja :: 1,1 (could be used for selecting subblock of a)
657 : !> \param alpha ...
658 : !> \param matrix_a ...
659 : !> \param ia ...
660 : !> \param ja ...
661 : !> \param beta ...
662 : !> \param matrix_c ...
663 : !> \author Matthias Krack
664 : ! **************************************************************************************************
665 11778 : SUBROUTINE cp_fm_syrk(uplo, trans, k, alpha, matrix_a, ia, ja, beta, matrix_c)
666 : CHARACTER(LEN=1), INTENT(IN) :: uplo, trans
667 : INTEGER, INTENT(IN) :: k
668 : REAL(KIND=dp), INTENT(IN) :: alpha
669 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a
670 : INTEGER, INTENT(IN) :: ia, ja
671 : REAL(KIND=dp), INTENT(IN) :: beta
672 : TYPE(cp_fm_type), INTENT(IN) :: matrix_c
673 :
674 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_syrk'
675 :
676 : INTEGER :: handle, n
677 11778 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, c
678 : #if defined(__parallel)
679 : INTEGER, DIMENSION(9) :: desca, descc
680 : #else
681 : INTEGER :: lda, ldc
682 : #endif
683 : #if defined (__HAS_IEEE_EXCEPTIONS)
684 : LOGICAL, DIMENSION(5) :: halt
685 : #endif
686 :
687 11778 : CALL timeset(routineN, handle)
688 :
689 11778 : n = matrix_c%matrix_struct%nrow_global
690 :
691 11778 : a => matrix_a%local_data
692 11778 : c => matrix_c%local_data
693 :
694 : #if defined (__HAS_IEEE_EXCEPTIONS)
695 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
696 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
697 : #endif
698 : #if defined(__parallel)
699 117780 : desca(:) = matrix_a%matrix_struct%descriptor(:)
700 117780 : descc(:) = matrix_c%matrix_struct%descriptor(:)
701 :
702 11778 : CALL pdsyrk(uplo, trans, n, k, alpha, a(1, 1), ia, ja, desca, beta, c(1, 1), 1, 1, descc)
703 : #else
704 : lda = SIZE(a, 1)
705 : ldc = SIZE(c, 1)
706 :
707 : CALL dsyrk(uplo, trans, n, k, alpha, a(ia, ja), lda, beta, c(1, 1), ldc)
708 : #endif
709 : #if defined (__HAS_IEEE_EXCEPTIONS)
710 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
711 : #endif
712 11778 : CALL timestop(handle)
713 :
714 11778 : END SUBROUTINE cp_fm_syrk
715 :
716 : ! **************************************************************************************************
717 : !> \brief computes the schur product of two matrices
718 : !> c_ij = a_ij * b_ij
719 : !> \param matrix_a ...
720 : !> \param matrix_b ...
721 : !> \param matrix_c ...
722 : !> \author Joost VandeVondele
723 : ! **************************************************************************************************
724 13066 : SUBROUTINE cp_fm_schur_product(matrix_a, matrix_b, matrix_c)
725 :
726 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a, matrix_b, matrix_c
727 :
728 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_schur_product'
729 :
730 : INTEGER :: handle, icol_local, irow_local, mypcol, &
731 : myprow, ncol_local, &
732 : nrow_local
733 13066 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, c
734 : TYPE(cp_blacs_env_type), POINTER :: context
735 :
736 13066 : CALL timeset(routineN, handle)
737 :
738 13066 : context => matrix_a%matrix_struct%context
739 13066 : myprow = context%mepos(1)
740 13066 : mypcol = context%mepos(2)
741 :
742 13066 : a => matrix_a%local_data
743 13066 : b => matrix_b%local_data
744 13066 : c => matrix_c%local_data
745 :
746 13066 : nrow_local = matrix_a%matrix_struct%nrow_locals(myprow)
747 13066 : ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol)
748 :
749 117646 : DO icol_local = 1, ncol_local
750 6743576 : DO irow_local = 1, nrow_local
751 6730510 : c(irow_local, icol_local) = a(irow_local, icol_local)*b(irow_local, icol_local)
752 : END DO
753 : END DO
754 :
755 13066 : CALL timestop(handle)
756 :
757 13066 : END SUBROUTINE cp_fm_schur_product
758 :
759 : ! **************************************************************************************************
760 : !> \brief returns the trace of matrix_a^T matrix_b, i.e
761 : !> sum_{i,j}(matrix_a(i,j)*matrix_b(i,j))
762 : !> \param matrix_a a matrix
763 : !> \param matrix_b another matrix
764 : !> \param trace ...
765 : !> \par History
766 : !> 11.06.2001 Creation (Matthias Krack)
767 : !> 12.2002 added doc [fawzi]
768 : !> \author Matthias Krack
769 : !> \note
770 : !> note the transposition of matrix_a!
771 : ! **************************************************************************************************
772 286702 : SUBROUTINE cp_fm_trace_a0b0t0(matrix_a, matrix_b, trace)
773 :
774 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a, matrix_b
775 : REAL(KIND=dp), INTENT(OUT) :: trace
776 :
777 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_trace_a0b0t0'
778 :
779 : INTEGER :: handle, mypcol, myprow, ncol_local, &
780 : nrow_local
781 286702 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
782 : TYPE(cp_blacs_env_type), POINTER :: context
783 : TYPE(mp_comm_type) :: group
784 :
785 286702 : CALL timeset(routineN, handle)
786 :
787 286702 : context => matrix_a%matrix_struct%context
788 286702 : myprow = context%mepos(1)
789 286702 : mypcol = context%mepos(2)
790 :
791 286702 : group = matrix_a%matrix_struct%para_env
792 :
793 286702 : a => matrix_a%local_data
794 286702 : b => matrix_b%local_data
795 :
796 286702 : nrow_local = MIN(matrix_a%matrix_struct%nrow_locals(myprow), matrix_b%matrix_struct%nrow_locals(myprow))
797 286702 : ncol_local = MIN(matrix_a%matrix_struct%ncol_locals(mypcol), matrix_b%matrix_struct%ncol_locals(mypcol))
798 :
799 : ! cries for an accurate_dot_product
800 : trace = accurate_dot_product(a(1:nrow_local, 1:ncol_local), &
801 286702 : b(1:nrow_local, 1:ncol_local))
802 :
803 286702 : CALL group%sum(trace)
804 :
805 286702 : CALL timestop(handle)
806 :
807 286702 : END SUBROUTINE cp_fm_trace_a0b0t0
808 :
809 : #:mute
810 : #:set types = [("cp_fm_type", "a", ""), ("cp_fm_p_type", "p","%matrix")]
811 : #:endmute
812 :
813 : ! **************************************************************************************************
814 : !> \brief Compute trace(k) = Tr (matrix_a(k)^T matrix_b) for each pair of matrices A_k and B.
815 : !> \param matrix_a list of A matrices
816 : !> \param matrix_b B matrix
817 : !> \param trace computed traces
818 : !> \par History
819 : !> * 08.2018 forked from cp_fm_trace() [Sergey Chulkov]
820 : !> \note \parblock
821 : !> Computing the trace requires collective communication between involved MPI processes
822 : !> that implies a synchronisation point between them. The aim of this subroutine is to reduce
823 : !> the amount of time wasted in such synchronisation by performing one large collective
824 : !> operation which involves all the matrices in question.
825 : !>
826 : !> The subroutine's suffix reflects dimensionality of dummy arrays; 'a1b0t1' means that
827 : !> the dummy variables 'matrix_a' and 'trace' are 1-dimensional arrays, while the variable
828 : !> 'matrix_b' is a single matrix.
829 : !> \endparblock
830 : ! **************************************************************************************************
831 : #:for longname, shortname, appendix in types
832 92473 : SUBROUTINE cp_fm_trace_a1b0t1_${shortname}$ (matrix_a, matrix_b, trace)
833 : TYPE(${longname}$), DIMENSION(:), INTENT(IN) :: matrix_a
834 : TYPE(cp_fm_type), INTENT(IN) :: matrix_b
835 : REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: trace
836 :
837 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_trace_a1b0t1_${shortname}$'
838 :
839 : INTEGER :: handle, imatrix, n_matrices, &
840 : ncols_local, nrows_local
841 92473 : REAL(kind=dp), DIMENSION(:, :), POINTER :: ldata_a, ldata_b
842 : TYPE(mp_comm_type) :: group
843 :
844 92473 : CALL timeset(routineN, handle)
845 :
846 92473 : n_matrices = SIZE(trace)
847 92473 : CPASSERT(SIZE(matrix_a) == n_matrices)
848 :
849 92473 : CALL cp_fm_get_info(matrix_b, nrow_local=nrows_local, ncol_local=ncols_local)
850 :
851 92473 : ldata_b => matrix_b%local_data(1:nrows_local, 1:ncols_local)
852 :
853 : !$OMP PARALLEL DO DEFAULT(NONE), &
854 : !$OMP PRIVATE(imatrix, ldata_a), &
855 : !$OMP SHARED(ldata_b, matrix_a, matrix_b), &
856 92473 : !$OMP SHARED(ncols_local, nrows_local, n_matrices, trace)
857 :
858 : DO imatrix = 1, n_matrices
859 :
860 : ! assume that the matrices A(i) and B have identical shapes and distribution schemes
861 : ldata_a => matrix_a(imatrix) ${appendix}$%local_data(1:nrows_local, 1:ncols_local)
862 : trace(imatrix) = accurate_dot_product(ldata_a, ldata_b)
863 : END DO
864 : !$OMP END PARALLEL DO
865 :
866 92473 : group = matrix_b%matrix_struct%para_env
867 743077 : CALL group%sum(trace)
868 :
869 92473 : CALL timestop(handle)
870 92473 : END SUBROUTINE cp_fm_trace_a1b0t1_${shortname}$
871 : #:endfor
872 :
873 : ! **************************************************************************************************
874 : !> \brief Compute trace(k) = Tr (matrix_a(k)^T matrix_b(k)) for each pair of matrices A_k and B_k.
875 : !> \param matrix_a list of A matrices
876 : !> \param matrix_b list of B matrices
877 : !> \param trace computed traces
878 : !> \param accurate ...
879 : !> \par History
880 : !> * 11.2016 forked from cp_fm_trace() [Sergey Chulkov]
881 : !> \note \parblock
882 : !> Computing the trace requires collective communication between involved MPI processes
883 : !> that implies a synchronisation point between them. The aim of this subroutine is to reduce
884 : !> the amount of time wasted in such synchronisation by performing one large collective
885 : !> operation which involves all the matrices in question.
886 : !>
887 : !> The subroutine's suffix reflects dimensionality of dummy arrays; 'a1b1t1' means that
888 : !> all dummy variables (matrix_a, matrix_b, and trace) are 1-dimensional arrays.
889 : !> \endparblock
890 : ! **************************************************************************************************
891 : #:for longname1, shortname1, appendix1 in types
892 : #:for longname2, shortname2, appendix2 in types
893 187268 : SUBROUTINE cp_fm_trace_a1b1t1_${shortname1}$${shortname2}$ (matrix_a, matrix_b, trace, accurate)
894 : TYPE(${longname1}$), DIMENSION(:), INTENT(IN) :: matrix_a
895 : TYPE(${longname2}$), DIMENSION(:), INTENT(IN) :: matrix_b
896 : REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: trace
897 : LOGICAL, INTENT(IN), OPTIONAL :: accurate
898 :
899 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_trace_a1b1t1_${shortname1}$${shortname2}$'
900 :
901 : INTEGER :: handle, imatrix, n_matrices, &
902 : ncols_local, nrows_local
903 : LOGICAL :: use_accurate_sum
904 187268 : REAL(kind=dp), DIMENSION(:, :), POINTER :: ldata_a, ldata_b
905 : TYPE(mp_comm_type) :: group
906 :
907 187268 : CALL timeset(routineN, handle)
908 :
909 187268 : n_matrices = SIZE(trace)
910 187268 : CPASSERT(SIZE(matrix_a) == n_matrices)
911 187268 : CPASSERT(SIZE(matrix_b) == n_matrices)
912 :
913 187268 : use_accurate_sum = .TRUE.
914 187268 : IF (PRESENT(accurate)) use_accurate_sum = accurate
915 :
916 : !$OMP PARALLEL DO DEFAULT(NONE), &
917 : !$OMP PRIVATE(imatrix, ldata_a, ldata_b, ncols_local), &
918 : !$OMP PRIVATE(nrows_local), &
919 187268 : !$OMP SHARED(matrix_a, matrix_b, n_matrices, trace, use_accurate_sum)
920 : DO imatrix = 1, n_matrices
921 : CALL cp_fm_get_info(matrix_a(imatrix) ${appendix1}$, nrow_local=nrows_local, ncol_local=ncols_local)
922 :
923 : ! assume that the matrices A(i) and B(i) have identical shapes and distribution schemes
924 : ldata_a => matrix_a(imatrix) ${appendix1}$%local_data(1:nrows_local, 1:ncols_local)
925 : ldata_b => matrix_b(imatrix) ${appendix2}$%local_data(1:nrows_local, 1:ncols_local)
926 : IF (use_accurate_sum) THEN
927 : trace(imatrix) = accurate_dot_product(ldata_a, ldata_b)
928 : ELSE
929 : trace(imatrix) = SUM(ldata_a*ldata_b)
930 : END IF
931 : END DO
932 : !$OMP END PARALLEL DO
933 :
934 187268 : group = matrix_a(1) ${appendix1}$%matrix_struct%para_env
935 650832 : CALL group%sum(trace)
936 :
937 187268 : CALL timestop(handle)
938 187268 : END SUBROUTINE cp_fm_trace_a1b1t1_${shortname1}$${shortname2}$
939 : #:endfor
940 : #:endfor
941 :
942 : ! **************************************************************************************************
943 : !> \brief Compute trace(i,j) = \sum_k Tr (matrix_a(k,i)^T matrix_b(k,j)).
944 : !> \param matrix_a list of A matrices
945 : !> \param matrix_b list of B matrices
946 : !> \param trace computed traces
947 : !> \param accurate ...
948 : ! **************************************************************************************************
949 : #:for longname1, shortname1, appendix1 in types
950 : #:for longname2, shortname2, appendix2 in types
951 17584 : SUBROUTINE cp_fm_contracted_trace_a2b2t2_${shortname1}$${shortname2}$ (matrix_a, matrix_b, trace, accurate)
952 : TYPE(${longname1}$), DIMENSION(:, :), INTENT(IN) :: matrix_a
953 : TYPE(${longname2}$), DIMENSION(:, :), INTENT(IN) :: matrix_b
954 : REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: trace
955 : LOGICAL, INTENT(IN), OPTIONAL :: accurate
956 :
957 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_contracted_trace_a2b2t2_${shortname1}$${shortname2}$'
958 :
959 : INTEGER :: handle, ia, ib, iz, na, nb, ncols_local, &
960 : nrows_local, nz
961 : INTEGER(kind=int_8) :: ib8, itrace, na8, ntraces
962 : LOGICAL :: use_accurate_sum
963 : REAL(kind=dp) :: t
964 17584 : REAL(kind=dp), DIMENSION(:, :), POINTER :: ldata_a, ldata_b
965 : TYPE(mp_comm_type) :: group
966 :
967 17584 : CALL timeset(routineN, handle)
968 :
969 17584 : nz = SIZE(matrix_a, 1)
970 17584 : CPASSERT(SIZE(matrix_b, 1) == nz)
971 :
972 17584 : na = SIZE(matrix_a, 2)
973 17584 : nb = SIZE(matrix_b, 2)
974 17584 : CPASSERT(SIZE(trace, 1) == na)
975 17584 : CPASSERT(SIZE(trace, 2) == nb)
976 :
977 17584 : use_accurate_sum = .TRUE.
978 17584 : IF (PRESENT(accurate)) use_accurate_sum = accurate
979 :
980 : ! here we use one running index (itrace) instead of two (ia, ib) in order to
981 : ! improve load balance between shared-memory threads
982 17584 : ntraces = na*nb
983 17584 : na8 = INT(na, kind=int_8)
984 :
985 : !$OMP PARALLEL DO DEFAULT(NONE), &
986 : !$OMP PRIVATE(ia, ib, ib8, itrace, iz, ldata_a, ldata_b, ncols_local), &
987 : !$OMP PRIVATE(nrows_local, t), &
988 17584 : !$OMP SHARED(matrix_a, matrix_b, na, na8, nb, ntraces, nz, trace, use_accurate_sum)
989 : DO itrace = 1, ntraces
990 : ib8 = (itrace - 1)/na8
991 : ia = INT(itrace - ib8*na8)
992 : ib = INT(ib8) + 1
993 :
994 : t = 0.0_dp
995 : DO iz = 1, nz
996 : CALL cp_fm_get_info(matrix_a(iz, ia) ${appendix1}$, nrow_local=nrows_local, ncol_local=ncols_local)
997 :
998 : ! assume that the matrices A(iz, ia) and B(iz, ib) have identical shapes and distribution schemes
999 : ldata_a => matrix_a(iz, ia) ${appendix1}$%local_data(1:nrows_local, 1:ncols_local)
1000 : ldata_b => matrix_b(iz, ib) ${appendix2}$%local_data(1:nrows_local, 1:ncols_local)
1001 : IF (use_accurate_sum) THEN
1002 : t = t + accurate_dot_product(ldata_a, ldata_b)
1003 : ELSE
1004 : t = t + SUM(ldata_a*ldata_b)
1005 : END IF
1006 : END DO
1007 : trace(ia, ib) = t
1008 : END DO
1009 : !$OMP END PARALLEL DO
1010 :
1011 17584 : group = matrix_a(1, 1) ${appendix1}$%matrix_struct%para_env
1012 833876 : CALL group%sum(trace)
1013 :
1014 17584 : CALL timestop(handle)
1015 17584 : END SUBROUTINE cp_fm_contracted_trace_a2b2t2_${shortname1}$${shortname2}$
1016 : #:endfor
1017 : #:endfor
1018 :
1019 : ! **************************************************************************************************
1020 : !> \brief multiplies in place by a triangular matrix:
1021 : !> matrix_b = alpha op(triangular_matrix) matrix_b
1022 : !> or (if side='R')
1023 : !> matrix_b = alpha matrix_b op(triangular_matrix)
1024 : !> op(triangular_matrix) is:
1025 : !> triangular_matrix (if transpose_tr=.false. and invert_tr=.false.)
1026 : !> triangular_matrix^T (if transpose_tr=.true. and invert_tr=.false.)
1027 : !> triangular_matrix^(-1) (if transpose_tr=.false. and invert_tr=.true.)
1028 : !> triangular_matrix^(-T) (if transpose_tr=.true. and invert_tr=.true.)
1029 : !> \param triangular_matrix the triangular matrix that multiplies the other
1030 : !> \param matrix_b the matrix that gets multiplied and stores the result
1031 : !> \param side on which side of matrix_b stays op(triangular_matrix)
1032 : !> (defaults to 'L')
1033 : !> \param transpose_tr if the triangular matrix should be transposed
1034 : !> (defaults to false)
1035 : !> \param invert_tr if the triangular matrix should be inverted
1036 : !> (defaults to false)
1037 : !> \param uplo_tr if triangular_matrix is stored in the upper ('U') or
1038 : !> lower ('L') triangle (defaults to 'U')
1039 : !> \param unit_diag_tr if the diagonal elements of triangular_matrix should
1040 : !> be assumed to be 1 (defaults to false)
1041 : !> \param n_rows the number of rows of the result (defaults to
1042 : !> size(matrix_b,1))
1043 : !> \param n_cols the number of columns of the result (defaults to
1044 : !> size(matrix_b,2))
1045 : !> \param alpha ...
1046 : !> \par History
1047 : !> 08.2002 created [fawzi]
1048 : !> \author Fawzi Mohamed
1049 : !> \note
1050 : !> needs an mpi env
1051 : ! **************************************************************************************************
1052 141232 : SUBROUTINE cp_fm_triangular_multiply(triangular_matrix, matrix_b, side, &
1053 : transpose_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, &
1054 : alpha)
1055 : TYPE(cp_fm_type), INTENT(IN) :: triangular_matrix, matrix_b
1056 : CHARACTER, INTENT(IN), OPTIONAL :: side
1057 : LOGICAL, INTENT(IN), OPTIONAL :: transpose_tr, invert_tr
1058 : CHARACTER, INTENT(IN), OPTIONAL :: uplo_tr
1059 : LOGICAL, INTENT(IN), OPTIONAL :: unit_diag_tr
1060 : INTEGER, INTENT(IN), OPTIONAL :: n_rows, n_cols
1061 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: alpha
1062 :
1063 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_triangular_multiply'
1064 :
1065 : CHARACTER :: side_char, transa, unit_diag, uplo
1066 : INTEGER :: handle, mdim, m, n
1067 : LOGICAL :: invert
1068 : REAL(KIND=dp) :: al
1069 :
1070 70616 : CALL timeset(routineN, handle)
1071 70616 : side_char = 'L'
1072 70616 : unit_diag = 'N'
1073 70616 : uplo = 'U'
1074 70616 : transa = 'N'
1075 70616 : invert = .FALSE.
1076 70616 : al = 1.0_dp
1077 70616 : CALL cp_fm_get_info(matrix_b, nrow_global=m, ncol_global=n)
1078 70616 : IF (PRESENT(side)) side_char = side
1079 : mdim = MERGE(1, 2, 'L' == side_char)
1080 70616 : IF (PRESENT(invert_tr)) invert = invert_tr
1081 70616 : IF (PRESENT(uplo_tr)) uplo = uplo_tr
1082 70616 : IF (PRESENT(unit_diag_tr)) THEN
1083 0 : IF (unit_diag_tr) THEN
1084 0 : unit_diag = 'U'
1085 : ELSE
1086 : unit_diag = 'N'
1087 : END IF
1088 : END IF
1089 70616 : IF (PRESENT(transpose_tr)) THEN
1090 6666 : IF (transpose_tr) THEN
1091 3028 : transa = 'T'
1092 : ELSE
1093 : transa = 'N'
1094 : END IF
1095 : END IF
1096 70616 : IF (PRESENT(alpha)) al = alpha
1097 70616 : IF (PRESENT(n_rows)) m = n_rows
1098 70616 : IF (PRESENT(n_cols)) n = n_cols
1099 :
1100 70616 : IF (invert) THEN
1101 :
1102 : #if defined(__parallel)
1103 : CALL pdtrsm(side_char, uplo, transa, unit_diag, m, n, al, &
1104 : triangular_matrix%local_data(1, 1), 1, 1, &
1105 : triangular_matrix%matrix_struct%descriptor, &
1106 : matrix_b%local_data(1, 1), 1, 1, &
1107 53150 : matrix_b%matrix_struct%descriptor(1))
1108 : #else
1109 : CALL dtrsm(side_char, uplo, transa, unit_diag, m, n, al, &
1110 : triangular_matrix%local_data(1, 1), &
1111 : SIZE(triangular_matrix%local_data, mdim), &
1112 : matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1))
1113 : #endif
1114 :
1115 : ELSE
1116 :
1117 : #if defined(__parallel)
1118 : CALL pdtrmm(side_char, uplo, transa, unit_diag, m, n, al, &
1119 : triangular_matrix%local_data(1, 1), 1, 1, &
1120 : triangular_matrix%matrix_struct%descriptor, &
1121 : matrix_b%local_data(1, 1), 1, 1, &
1122 17466 : matrix_b%matrix_struct%descriptor(1))
1123 : #else
1124 : CALL dtrmm(side_char, uplo, transa, unit_diag, m, n, al, &
1125 : triangular_matrix%local_data(1, 1), &
1126 : SIZE(triangular_matrix%local_data, mdim), &
1127 : matrix_b%local_data(1, 1), SIZE(matrix_b%local_data, 1))
1128 : #endif
1129 :
1130 : END IF
1131 :
1132 70616 : CALL timestop(handle)
1133 70616 : END SUBROUTINE cp_fm_triangular_multiply
1134 :
1135 : ! **************************************************************************************************
1136 : !> \brief scales a matrix
1137 : !> matrix_a = alpha * matrix_b
1138 : !> \param alpha ...
1139 : !> \param matrix_a ...
1140 : !> \note
1141 : !> use cp_fm_set_all to zero (avoids problems with nan)
1142 : ! **************************************************************************************************
1143 231780 : SUBROUTINE cp_fm_scale(alpha, matrix_a)
1144 : REAL(KIND=dp), INTENT(IN) :: alpha
1145 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a
1146 :
1147 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_scale'
1148 :
1149 : INTEGER :: handle, size_a
1150 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
1151 :
1152 231780 : CALL timeset(routineN, handle)
1153 :
1154 : NULLIFY (a)
1155 :
1156 231780 : a => matrix_a%local_data
1157 231780 : size_a = SIZE(a, 1)*SIZE(a, 2)
1158 :
1159 231780 : CALL DSCAL(size_a, alpha, a, 1)
1160 :
1161 231780 : CALL timestop(handle)
1162 :
1163 231780 : END SUBROUTINE cp_fm_scale
1164 :
1165 : ! **************************************************************************************************
1166 : !> \brief transposes a matrix
1167 : !> matrixt = matrix ^ T
1168 : !> \param matrix ...
1169 : !> \param matrixt ...
1170 : !> \note
1171 : !> all matrix elements are transposed (see cp_fm_uplo_to_full to symmetrise a matrix)
1172 : ! **************************************************************************************************
1173 40688 : SUBROUTINE cp_fm_transpose(matrix, matrixt)
1174 : TYPE(cp_fm_type), INTENT(IN) :: matrix, matrixt
1175 :
1176 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_transpose'
1177 :
1178 : INTEGER :: handle, ncol_global, &
1179 : nrow_global, ncol_globalt, nrow_globalt
1180 20344 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, c
1181 : #if defined(__parallel)
1182 : INTEGER, DIMENSION(9) :: desca, descc
1183 : #elif !defined(__MKL)
1184 : INTEGER :: i, j
1185 : #endif
1186 :
1187 20344 : nrow_global = matrix%matrix_struct%nrow_global
1188 20344 : ncol_global = matrix%matrix_struct%ncol_global
1189 20344 : nrow_globalt = matrixt%matrix_struct%nrow_global
1190 20344 : ncol_globalt = matrixt%matrix_struct%ncol_global
1191 0 : CPASSERT(nrow_global == ncol_globalt)
1192 20344 : CPASSERT(nrow_globalt == ncol_global)
1193 :
1194 20344 : CALL timeset(routineN, handle)
1195 :
1196 20344 : a => matrix%local_data
1197 20344 : c => matrixt%local_data
1198 :
1199 : #if defined(__parallel)
1200 203440 : desca(:) = matrix%matrix_struct%descriptor(:)
1201 203440 : descc(:) = matrixt%matrix_struct%descriptor(:)
1202 20344 : CALL pdtran(ncol_global, nrow_global, 1.0_dp, a(1, 1), 1, 1, desca, 0.0_dp, c(1, 1), 1, 1, descc)
1203 : #elif defined(__MKL)
1204 : CALL mkl_domatcopy('C', 'T', nrow_global, ncol_global, 1.0_dp, a(1, 1), nrow_global, c(1, 1), ncol_global)
1205 : #else
1206 : DO j = 1, ncol_global
1207 : DO i = 1, nrow_global
1208 : c(j, i) = a(i, j)
1209 : END DO
1210 : END DO
1211 : #endif
1212 20344 : CALL timestop(handle)
1213 :
1214 20344 : END SUBROUTINE cp_fm_transpose
1215 :
1216 : ! **************************************************************************************************
1217 : !> \brief given a triangular matrix according to uplo, computes the corresponding full matrix
1218 : !> \param matrix the triangular matrix as input, the full matrix as output
1219 : !> \param work a matrix of the same size as matrix
1220 : !> \param uplo triangular format; defaults to 'U'
1221 : !> \author Matthias Krack
1222 : !> \note
1223 : !> the opposite triangular part is irrelevant
1224 : ! **************************************************************************************************
1225 387746 : SUBROUTINE cp_fm_uplo_to_full(matrix, work, uplo)
1226 :
1227 : TYPE(cp_fm_type), INTENT(IN) :: matrix, work
1228 : CHARACTER, INTENT(IN), OPTIONAL :: uplo
1229 :
1230 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_uplo_to_full'
1231 :
1232 : CHARACTER :: myuplo
1233 : INTEGER :: handle, icol_global, irow_global, &
1234 : mypcol, myprow, ncol_global, &
1235 : nrow_global
1236 193873 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
1237 : TYPE(cp_blacs_env_type), POINTER :: context
1238 :
1239 : #if defined(__parallel)
1240 : INTEGER :: icol_local, irow_local, &
1241 : ncol_block, ncol_local, &
1242 : nrow_block, nrow_local
1243 : INTEGER, DIMENSION(9) :: desca, descc
1244 193873 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: c
1245 : #endif
1246 :
1247 193873 : myuplo = 'U'
1248 0 : IF (PRESENT(uplo)) myuplo = uplo
1249 :
1250 193873 : nrow_global = matrix%matrix_struct%nrow_global
1251 193873 : ncol_global = matrix%matrix_struct%ncol_global
1252 193873 : CPASSERT(nrow_global == ncol_global)
1253 193873 : nrow_global = work%matrix_struct%nrow_global
1254 193873 : ncol_global = work%matrix_struct%ncol_global
1255 193873 : CPASSERT(nrow_global == ncol_global)
1256 :
1257 193873 : CALL timeset(routineN, handle)
1258 :
1259 193873 : context => matrix%matrix_struct%context
1260 193873 : myprow = context%mepos(1)
1261 193873 : mypcol = context%mepos(2)
1262 :
1263 : #if defined(__parallel)
1264 :
1265 193873 : nrow_block = matrix%matrix_struct%nrow_block
1266 193873 : ncol_block = matrix%matrix_struct%ncol_block
1267 :
1268 193873 : nrow_local = matrix%matrix_struct%nrow_locals(myprow)
1269 193873 : ncol_local = matrix%matrix_struct%ncol_locals(mypcol)
1270 :
1271 193873 : a => work%local_data
1272 1938730 : desca(:) = work%matrix_struct%descriptor(:)
1273 193873 : c => matrix%local_data
1274 1938730 : descc(:) = matrix%matrix_struct%descriptor(:)
1275 :
1276 6026996 : DO icol_local = 1, ncol_local
1277 5833123 : icol_global = matrix%matrix_struct%col_indices(icol_local)
1278 245655146 : DO irow_local = 1, nrow_local
1279 239628150 : irow_global = matrix%matrix_struct%row_indices(irow_local)
1280 245461273 : IF (MERGE(irow_global > icol_global, irow_global < icol_global, (myuplo == "U") .OR. (myuplo == "u"))) THEN
1281 118230619 : c(irow_local, icol_local) = 0.0_dp
1282 121397531 : ELSE IF (irow_global == icol_global) THEN
1283 3166884 : c(irow_local, icol_local) = 0.5_dp*c(irow_local, icol_local)
1284 : END IF
1285 : END DO
1286 : END DO
1287 :
1288 6026996 : DO icol_local = 1, ncol_local
1289 245655146 : DO irow_local = 1, nrow_local
1290 245461273 : a(irow_local, icol_local) = c(irow_local, icol_local)
1291 : END DO
1292 : END DO
1293 :
1294 193873 : CALL pdtran(nrow_global, ncol_global, 1.0_dp, a(1, 1), 1, 1, desca, 1.0_dp, c(1, 1), 1, 1, descc)
1295 :
1296 : #else
1297 :
1298 : a => matrix%local_data
1299 :
1300 : IF ((myuplo == "U") .OR. (myuplo == "u")) THEN
1301 : DO irow_global = 1, nrow_global
1302 : DO icol_global = irow_global + 1, ncol_global
1303 : a(icol_global, irow_global) = a(irow_global, icol_global)
1304 : END DO
1305 : END DO
1306 : ELSE
1307 : DO icol_global = 1, ncol_global
1308 : DO irow_global = icol_global + 1, nrow_global
1309 : a(irow_global, icol_global) = a(icol_global, irow_global)
1310 : END DO
1311 : END DO
1312 : END IF
1313 :
1314 : #endif
1315 193873 : CALL timestop(handle)
1316 :
1317 193873 : END SUBROUTINE cp_fm_uplo_to_full
1318 :
1319 : ! **************************************************************************************************
1320 : !> \brief scales column i of matrix a with scaling(i)
1321 : !> \param matrixa ...
1322 : !> \param scaling : an array used for scaling the columns,
1323 : !> SIZE(scaling) determines the number of columns to be scaled
1324 : !> \author Joost VandeVondele
1325 : !> \note
1326 : !> this is very useful as a first step in the computation of C = sum_i alpha_i A_i transpose (A_i)
1327 : !> that is a rank-k update (cp_fm_syrk , cp_sm_plus_fm_fm_t)
1328 : !> this procedure can be up to 20 times faster than calling cp_fm_syrk n times
1329 : !> where every vector has a different prefactor
1330 : ! **************************************************************************************************
1331 122474 : SUBROUTINE cp_fm_column_scale(matrixa, scaling)
1332 : TYPE(cp_fm_type), INTENT(IN) :: matrixa
1333 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: scaling
1334 :
1335 : INTEGER :: k, mypcol, myprow, n, ncol_global, &
1336 : npcol, nprow
1337 122474 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
1338 : #if defined(__parallel)
1339 : INTEGER :: icol_global, icol_local, &
1340 : ipcol, iprow, irow_local
1341 : #else
1342 : INTEGER :: i
1343 : #endif
1344 :
1345 122474 : myprow = matrixa%matrix_struct%context%mepos(1)
1346 122474 : mypcol = matrixa%matrix_struct%context%mepos(2)
1347 122474 : nprow = matrixa%matrix_struct%context%num_pe(1)
1348 122474 : npcol = matrixa%matrix_struct%context%num_pe(2)
1349 :
1350 122474 : ncol_global = matrixa%matrix_struct%ncol_global
1351 :
1352 122474 : a => matrixa%local_data
1353 122474 : n = SIZE(a, 1)
1354 122474 : k = MIN(SIZE(scaling), ncol_global)
1355 :
1356 : #if defined(__parallel)
1357 :
1358 1237438 : DO icol_global = 1, k
1359 : CALL infog2l(1, icol_global, matrixa%matrix_struct%descriptor, &
1360 : nprow, npcol, myprow, mypcol, &
1361 1114964 : irow_local, icol_local, iprow, ipcol)
1362 1237438 : IF ((ipcol == mypcol)) THEN
1363 1114964 : CALL DSCAL(n, scaling(icol_global), a(:, icol_local), 1)
1364 : END IF
1365 : END DO
1366 : #else
1367 : DO i = 1, k
1368 : CALL DSCAL(n, scaling(i), a(:, i), 1)
1369 : END DO
1370 : #endif
1371 122474 : END SUBROUTINE cp_fm_column_scale
1372 :
1373 : ! **************************************************************************************************
1374 : !> \brief scales row i of matrix a with scaling(i)
1375 : !> \param matrixa ...
1376 : !> \param scaling : an array used for scaling the columns,
1377 : !> \author JGH
1378 : !> \note
1379 : ! **************************************************************************************************
1380 8214 : SUBROUTINE cp_fm_row_scale(matrixa, scaling)
1381 : TYPE(cp_fm_type), INTENT(IN) :: matrixa
1382 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: scaling
1383 :
1384 : INTEGER :: n, m, nrow_global, nrow_local, ncol_local
1385 8214 : INTEGER, DIMENSION(:), POINTER :: row_indices
1386 8214 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
1387 : #if defined(__parallel)
1388 : INTEGER :: irow_global, icol, irow
1389 : #else
1390 : INTEGER :: j
1391 : #endif
1392 :
1393 : CALL cp_fm_get_info(matrixa, row_indices=row_indices, nrow_global=nrow_global, &
1394 8214 : nrow_local=nrow_local, ncol_local=ncol_local)
1395 8214 : CPASSERT(SIZE(scaling) == nrow_global)
1396 :
1397 8214 : a => matrixa%local_data
1398 8214 : n = SIZE(a, 1)
1399 8214 : m = SIZE(a, 2)
1400 :
1401 : #if defined(__parallel)
1402 88030 : DO icol = 1, ncol_local
1403 6522512 : DO irow = 1, nrow_local
1404 6434482 : irow_global = row_indices(irow)
1405 6514298 : a(irow, icol) = scaling(irow_global)*a(irow, icol)
1406 : END DO
1407 : END DO
1408 : #else
1409 : DO j = 1, m
1410 : a(1:n, j) = scaling(1:n)*a(1:n, j)
1411 : END DO
1412 : #endif
1413 8214 : END SUBROUTINE cp_fm_row_scale
1414 :
1415 : ! **************************************************************************************************
1416 : !> \brief Inverts a cp_fm_type matrix, optionally returning the determinant of the input matrix
1417 : !> \param matrix_a the matrix to invert
1418 : !> \param matrix_inverse the inverse of matrix_a
1419 : !> \param det_a the determinant of matrix_a
1420 : !> \param eps_svd optional parameter to active SVD based inversion, singular values below eps_svd
1421 : !> are screened
1422 : !> \param eigval optionally return matrix eigenvalues/singular values
1423 : !> \par History
1424 : !> note of Jan Wilhelm (12.2015)
1425 : !> - computation of determinant corrected
1426 : !> - determinant only computed if det_a is present
1427 : !> 12.2016 added option to use SVD instead of LU [Nico Holmberg]
1428 : !> - Use cp_fm_get diag instead of n times cp_fm_get_element (A. Bussy)
1429 : !> \author Florian Schiffmann(02.2007)
1430 : ! **************************************************************************************************
1431 1854 : SUBROUTINE cp_fm_invert(matrix_a, matrix_inverse, det_a, eps_svd, eigval)
1432 :
1433 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a, matrix_inverse
1434 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: det_a
1435 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_svd
1436 : REAL(KIND=dp), DIMENSION(:), POINTER, &
1437 : INTENT(INOUT), OPTIONAL :: eigval
1438 :
1439 : INTEGER :: n
1440 1854 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
1441 : REAL(KIND=dp) :: determinant, my_eps_svd
1442 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
1443 : TYPE(cp_fm_type) :: matrix_lu
1444 :
1445 : #if defined(__parallel)
1446 : TYPE(cp_fm_type) :: u, vt, sigma, inv_sigma_ut
1447 : TYPE(mp_comm_type) :: group
1448 : INTEGER :: i, info, liwork, lwork, exponent_of_minus_one
1449 : INTEGER, DIMENSION(9) :: desca
1450 : LOGICAL :: quenched
1451 : REAL(KIND=dp) :: alpha, beta
1452 1854 : REAL(KIND=dp), DIMENSION(:), POINTER :: diag
1453 1854 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
1454 : #else
1455 : LOGICAL :: sign
1456 : REAL(KIND=dp) :: eps1
1457 : #endif
1458 :
1459 1854 : my_eps_svd = 0.0_dp
1460 480 : IF (PRESENT(eps_svd)) my_eps_svd = eps_svd
1461 :
1462 : CALL cp_fm_create(matrix=matrix_lu, &
1463 : matrix_struct=matrix_a%matrix_struct, &
1464 1854 : name="A_lu"//TRIM(ADJUSTL(cp_to_string(1)))//"MATRIX")
1465 1854 : CALL cp_fm_to_fm(matrix_a, matrix_lu)
1466 :
1467 1854 : a => matrix_lu%local_data
1468 1854 : n = matrix_lu%matrix_struct%nrow_global
1469 5562 : ALLOCATE (ipivot(n + matrix_a%matrix_struct%nrow_block))
1470 1854 : ipivot(:) = 0
1471 : #if defined(__parallel)
1472 1854 : IF (my_eps_svd == 0.0_dp) THEN
1473 : ! Use LU decomposition
1474 : lwork = 3*n
1475 : liwork = 3*n
1476 18240 : desca(:) = matrix_lu%matrix_struct%descriptor(:)
1477 1824 : CALL pdgetrf(n, n, a, 1, 1, desca, ipivot, info)
1478 :
1479 1824 : IF (PRESENT(det_a) .OR. PRESENT(eigval)) THEN
1480 :
1481 1134 : ALLOCATE (diag(n))
1482 928 : diag(:) = 0.0_dp
1483 378 : CALL cp_fm_get_diag(matrix_lu, diag)
1484 :
1485 378 : exponent_of_minus_one = 0
1486 378 : determinant = 1.0_dp
1487 928 : DO i = 1, n
1488 550 : determinant = determinant*diag(i)
1489 928 : IF (ipivot(i) /= i) THEN
1490 224 : exponent_of_minus_one = exponent_of_minus_one + 1
1491 : END IF
1492 : END DO
1493 378 : IF (PRESENT(eigval)) THEN
1494 0 : CPASSERT(.NOT. ASSOCIATED(eigval))
1495 0 : ALLOCATE (eigval(n))
1496 0 : eigval(:) = diag
1497 : END IF
1498 378 : DEALLOCATE (diag)
1499 :
1500 378 : group = matrix_lu%matrix_struct%para_env
1501 378 : CALL group%sum(exponent_of_minus_one)
1502 :
1503 378 : determinant = determinant*(-1.0_dp)**exponent_of_minus_one
1504 :
1505 : END IF
1506 :
1507 1824 : alpha = 0.0_dp
1508 1824 : beta = 1.0_dp
1509 1824 : CALL cp_fm_set_all(matrix_inverse, alpha, beta)
1510 1824 : CALL pdgetrs('N', n, n, matrix_lu%local_data, 1, 1, desca, ipivot, matrix_inverse%local_data, 1, 1, desca, info)
1511 : ELSE
1512 : ! Use singular value decomposition
1513 : CALL cp_fm_create(matrix=u, &
1514 : matrix_struct=matrix_a%matrix_struct, &
1515 30 : name="LEFT_SINGULAR_MATRIX")
1516 30 : CALL cp_fm_set_all(u, alpha=0.0_dp)
1517 : CALL cp_fm_create(matrix=vt, &
1518 : matrix_struct=matrix_a%matrix_struct, &
1519 30 : name="RIGHT_SINGULAR_MATRIX")
1520 30 : CALL cp_fm_set_all(vt, alpha=0.0_dp)
1521 90 : ALLOCATE (diag(n))
1522 96 : diag(:) = 0.0_dp
1523 300 : desca(:) = matrix_lu%matrix_struct%descriptor(:)
1524 30 : ALLOCATE (work(1))
1525 : ! Workspace query
1526 30 : lwork = -1
1527 : CALL pdgesvd('V', 'V', n, n, matrix_lu%local_data, 1, 1, desca, diag, u%local_data, &
1528 30 : 1, 1, desca, vt%local_data, 1, 1, desca, work, lwork, info)
1529 30 : lwork = INT(work(1))
1530 30 : DEALLOCATE (work)
1531 90 : ALLOCATE (work(lwork))
1532 : ! SVD
1533 : CALL pdgesvd('V', 'V', n, n, matrix_lu%local_data, 1, 1, desca, diag, u%local_data, &
1534 30 : 1, 1, desca, vt%local_data, 1, 1, desca, work, lwork, info)
1535 : ! info == n+1 implies homogeneity error when the number of procs is large
1536 : ! this likely isnt a problem, but maybe we should handle it separately
1537 30 : IF (info /= 0 .AND. info /= n + 1) &
1538 0 : CPABORT("Singular value decomposition of matrix failed.")
1539 : ! (Pseudo)inverse and (pseudo)determinant
1540 : CALL cp_fm_create(matrix=sigma, &
1541 : matrix_struct=matrix_a%matrix_struct, &
1542 30 : name="SINGULAR_VALUE_MATRIX")
1543 30 : CALL cp_fm_set_all(sigma, alpha=0.0_dp)
1544 30 : determinant = 1.0_dp
1545 30 : quenched = .FALSE.
1546 30 : IF (PRESENT(eigval)) THEN
1547 28 : CPASSERT(.NOT. ASSOCIATED(eigval))
1548 84 : ALLOCATE (eigval(n))
1549 156 : eigval(:) = diag
1550 : END IF
1551 96 : DO i = 1, n
1552 66 : IF (diag(i) < my_eps_svd) THEN
1553 18 : diag(i) = 0.0_dp
1554 18 : quenched = .TRUE.
1555 : ELSE
1556 48 : determinant = determinant*diag(i)
1557 48 : diag(i) = 1.0_dp/diag(i)
1558 : END IF
1559 96 : CALL cp_fm_set_element(sigma, i, i, diag(i))
1560 : END DO
1561 30 : DEALLOCATE (diag)
1562 30 : IF (quenched) &
1563 : CALL cp_warn(__LOCATION__, &
1564 : "Linear dependencies were detected in the SVD inversion of matrix "//TRIM(ADJUSTL(matrix_a%name))// &
1565 12 : ". At least one singular value has been quenched.")
1566 : ! Sigma^-1 * U^T
1567 : CALL cp_fm_create(matrix=inv_sigma_ut, &
1568 : matrix_struct=matrix_a%matrix_struct, &
1569 30 : name="SINGULAR_VALUE_MATRIX")
1570 30 : CALL cp_fm_set_all(inv_sigma_ut, alpha=0.0_dp)
1571 : CALL pdgemm('N', 'T', n, n, n, 1.0_dp, sigma%local_data, 1, 1, desca, &
1572 30 : u%local_data, 1, 1, desca, 0.0_dp, inv_sigma_ut%local_data, 1, 1, desca)
1573 : ! A^-1 = V * (Sigma^-1 * U^T)
1574 30 : CALL cp_fm_set_all(matrix_inverse, alpha=0.0_dp)
1575 : CALL pdgemm('T', 'N', n, n, n, 1.0_dp, vt%local_data, 1, 1, desca, &
1576 30 : inv_sigma_ut%local_data, 1, 1, desca, 0.0_dp, matrix_inverse%local_data, 1, 1, desca)
1577 : ! Clean up
1578 30 : DEALLOCATE (work)
1579 30 : CALL cp_fm_release(u)
1580 30 : CALL cp_fm_release(vt)
1581 30 : CALL cp_fm_release(sigma)
1582 150 : CALL cp_fm_release(inv_sigma_ut)
1583 : END IF
1584 : #else
1585 : IF (my_eps_svd == 0.0_dp) THEN
1586 : sign = .TRUE.
1587 : CALL invert_matrix(matrix_a%local_data, matrix_inverse%local_data, &
1588 : eval_error=eps1)
1589 : CALL cp_fm_lu_decompose(matrix_lu, determinant, correct_sign=sign)
1590 : IF (PRESENT(eigval)) &
1591 : CALL cp_abort(__LOCATION__, &
1592 : "NYI. Eigenvalues not available for return without SCALAPACK.")
1593 : ELSE
1594 : CALL get_pseudo_inverse_svd(matrix_a%local_data, matrix_inverse%local_data, eps_svd, &
1595 : determinant, eigval)
1596 : END IF
1597 : #endif
1598 1854 : CALL cp_fm_release(matrix_lu)
1599 1854 : DEALLOCATE (ipivot)
1600 1854 : IF (PRESENT(det_a)) det_a = determinant
1601 1854 : END SUBROUTINE cp_fm_invert
1602 :
1603 : ! **************************************************************************************************
1604 : !> \brief inverts a triangular matrix
1605 : !> \param matrix_a ...
1606 : !> \param uplo_tr triangular format; defaults to 'U'
1607 : !> \author MI
1608 : ! **************************************************************************************************
1609 5880 : SUBROUTINE cp_fm_triangular_invert(matrix_a, uplo_tr)
1610 :
1611 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a
1612 : CHARACTER, INTENT(IN), OPTIONAL :: uplo_tr
1613 :
1614 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_triangular_invert'
1615 :
1616 : CHARACTER :: unit_diag, uplo
1617 : INTEGER :: handle, info, ncol_global
1618 5880 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
1619 : #if defined(__parallel)
1620 : INTEGER, DIMENSION(9) :: desca
1621 : #endif
1622 :
1623 5880 : CALL timeset(routineN, handle)
1624 :
1625 5880 : unit_diag = 'N'
1626 5880 : uplo = 'U'
1627 5880 : IF (PRESENT(uplo_tr)) uplo = uplo_tr
1628 :
1629 5880 : ncol_global = matrix_a%matrix_struct%ncol_global
1630 :
1631 5880 : a => matrix_a%local_data
1632 :
1633 : #if defined(__parallel)
1634 58800 : desca(:) = matrix_a%matrix_struct%descriptor(:)
1635 :
1636 5880 : CALL pdtrtri(uplo, unit_diag, ncol_global, a(1, 1), 1, 1, desca, info)
1637 :
1638 : #else
1639 : CALL dtrtri(uplo, unit_diag, ncol_global, a(1, 1), ncol_global, info)
1640 : #endif
1641 :
1642 5880 : CALL timestop(handle)
1643 5880 : END SUBROUTINE cp_fm_triangular_invert
1644 :
1645 : ! **************************************************************************************************
1646 : !> \brief performs a QR factorization of the input rectangular matrix A or of a submatrix of A
1647 : !> the computed triangular matrix R is in output of the submatrix sub(A) of size NxN
1648 : !> M and M give the dimension of the submatrix that has to be factorized (MxN) with M>N
1649 : !> \param matrix_a ...
1650 : !> \param matrix_r ...
1651 : !> \param nrow_fact ...
1652 : !> \param ncol_fact ...
1653 : !> \param first_row ...
1654 : !> \param first_col ...
1655 : !> \author MI
1656 : ! **************************************************************************************************
1657 23296 : SUBROUTINE cp_fm_qr_factorization(matrix_a, matrix_r, nrow_fact, ncol_fact, first_row, first_col, uplo)
1658 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a, matrix_r
1659 : INTEGER, INTENT(IN), OPTIONAL :: nrow_fact, ncol_fact, &
1660 : first_row, first_col
1661 : CHARACTER, INTENT(IN), OPTIONAL :: uplo
1662 :
1663 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_qr_factorization'
1664 :
1665 : CHARACTER :: myuplo
1666 : INTEGER :: handle, i, icol, info, irow, &
1667 : j, lda, lwork, ncol, &
1668 : ndim, nrow
1669 23296 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: tau, work
1670 23296 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: r_mat
1671 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
1672 : #if defined(__parallel)
1673 : INTEGER, DIMENSION(9) :: desca
1674 : #endif
1675 :
1676 23296 : CALL timeset(routineN, handle)
1677 :
1678 23296 : myuplo = 'U'
1679 23296 : IF (PRESENT(uplo)) myuplo = uplo
1680 :
1681 23296 : ncol = matrix_a%matrix_struct%ncol_global
1682 23296 : nrow = matrix_a%matrix_struct%nrow_global
1683 23296 : lda = nrow
1684 :
1685 23296 : a => matrix_a%local_data
1686 :
1687 23296 : IF (PRESENT(nrow_fact)) nrow = nrow_fact
1688 23296 : IF (PRESENT(ncol_fact)) ncol = ncol_fact
1689 23296 : irow = 1
1690 23296 : IF (PRESENT(first_row)) irow = first_row
1691 23296 : icol = 1
1692 23296 : IF (PRESENT(first_col)) icol = first_col
1693 :
1694 23296 : CPASSERT(nrow >= ncol)
1695 23296 : ndim = SIZE(a, 2)
1696 69888 : ALLOCATE (tau(ndim))
1697 :
1698 : #if defined(__parallel)
1699 :
1700 232960 : desca(:) = matrix_a%matrix_struct%descriptor(:)
1701 :
1702 23296 : lwork = -1
1703 69888 : ALLOCATE (work(2*ndim))
1704 23296 : CALL pdgeqrf(nrow, ncol, a, irow, icol, desca, tau, work, lwork, info)
1705 23296 : lwork = INT(work(1))
1706 23296 : DEALLOCATE (work)
1707 69888 : ALLOCATE (work(lwork))
1708 23296 : CALL pdgeqrf(nrow, ncol, a, irow, icol, desca, tau, work, lwork, info)
1709 :
1710 : #else
1711 : lwork = -1
1712 : ALLOCATE (work(2*ndim))
1713 : CALL dgeqrf(nrow, ncol, a, lda, tau, work, lwork, info)
1714 : lwork = INT(work(1))
1715 : DEALLOCATE (work)
1716 : ALLOCATE (work(lwork))
1717 : CALL dgeqrf(nrow, ncol, a, lda, tau, work, lwork, info)
1718 :
1719 : #endif
1720 :
1721 93184 : ALLOCATE (r_mat(ncol, ncol))
1722 23296 : CALL cp_fm_get_submatrix(matrix_a, r_mat, 1, 1, ncol, ncol)
1723 23296 : IF ((myuplo == "U") .OR. (myuplo == "u")) THEN
1724 56840 : DO i = 1, ncol
1725 76734 : DO j = i + 1, ncol
1726 53438 : r_mat(j, i) = 0.0_dp
1727 : END DO
1728 : END DO
1729 : ELSE
1730 0 : DO j = 1, ncol
1731 0 : DO i = j + 1, ncol
1732 0 : r_mat(i, j) = 0.0_dp
1733 : END DO
1734 : END DO
1735 : END IF
1736 23296 : CALL cp_fm_set_submatrix(matrix_r, r_mat, 1, 1, ncol, ncol)
1737 :
1738 23296 : DEALLOCATE (tau, work, r_mat)
1739 :
1740 23296 : CALL timestop(handle)
1741 :
1742 23296 : END SUBROUTINE cp_fm_qr_factorization
1743 :
1744 : ! **************************************************************************************************
1745 : !> \brief computes the the solution to A*b=A_general using lu decomposition
1746 : !> \param matrix_a input matrix; will be overwritten
1747 : !> \param general_a contains the result
1748 : !> \author Florian Schiffmann
1749 : ! **************************************************************************************************
1750 8094 : SUBROUTINE cp_fm_solve(matrix_a, general_a)
1751 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a, general_a
1752 :
1753 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_solve'
1754 :
1755 : INTEGER :: handle, info, n, nrhs
1756 8094 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot
1757 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, a_general
1758 : #if defined(__parallel)
1759 : INTEGER, DIMENSION(9) :: desca, descb
1760 : #else
1761 : INTEGER :: lda, ldb
1762 : #endif
1763 :
1764 8094 : CALL timeset(routineN, handle)
1765 :
1766 8094 : a => matrix_a%local_data
1767 8094 : a_general => general_a%local_data
1768 8094 : n = matrix_a%matrix_struct%nrow_global
1769 8094 : nrhs = general_a%matrix_struct%ncol_global
1770 24282 : ALLOCATE (ipivot(n + matrix_a%matrix_struct%nrow_block))
1771 :
1772 : #if defined(__parallel)
1773 80940 : desca(:) = matrix_a%matrix_struct%descriptor(:)
1774 80940 : descb(:) = general_a%matrix_struct%descriptor(:)
1775 8094 : CALL pdgetrf(n, n, a, 1, 1, desca, ipivot, info)
1776 : CALL pdgetrs("N", n, nrhs, a, 1, 1, desca, ipivot, a_general, &
1777 8094 : 1, 1, descb, info)
1778 :
1779 : #else
1780 : lda = SIZE(a, 1)
1781 : ldb = SIZE(a_general, 1)
1782 : CALL dgetrf(n, n, a, lda, ipivot, info)
1783 : CALL dgetrs("N", n, nrhs, a, lda, ipivot, a_general, ldb, info)
1784 :
1785 : #endif
1786 : ! info is allowed to be zero
1787 : ! this does just signal a zero diagonal element
1788 8094 : DEALLOCATE (ipivot)
1789 8094 : CALL timestop(handle)
1790 8094 : END SUBROUTINE
1791 :
1792 : ! **************************************************************************************************
1793 : !> \brief Convenience function. Computes the matrix multiplications needed
1794 : !> for the multiplication of complex matrices.
1795 : !> C = beta * C + alpha * ( A ** transa ) * ( B ** transb )
1796 : !> \param transa : 'N' -> normal 'T' -> transpose 'C' -> conjugate transpose
1797 : !> alpha,beta :: can be 0.0_dp and 1.0_dp
1798 : !> \param transb ...
1799 : !> \param m ...
1800 : !> \param n ...
1801 : !> \param k ...
1802 : !> \param alpha ...
1803 : !> \param A_re m x k matrix ( ! for transa = 'N'), real part
1804 : !> \param A_im m x k matrix ( ! for transa = 'N'), imaginary part
1805 : !> \param B_re k x n matrix ( ! for transa = 'N'), real part
1806 : !> \param B_im k x n matrix ( ! for transa = 'N'), imaginary part
1807 : !> \param beta ...
1808 : !> \param C_re m x n matrix, real part
1809 : !> \param C_im m x n matrix, imaginary part
1810 : !> \param a_first_col ...
1811 : !> \param a_first_row ...
1812 : !> \param b_first_col : the k x n matrix starts at col b_first_col of matrix_b (avoid usage)
1813 : !> \param b_first_row ...
1814 : !> \param c_first_col ...
1815 : !> \param c_first_row ...
1816 : !> \author Samuel Andermatt
1817 : !> \note
1818 : !> C should have no overlap with A, B
1819 : ! **************************************************************************************************
1820 4748 : SUBROUTINE cp_complex_fm_gemm(transa, transb, m, n, k, alpha, A_re, A_im, B_re, B_im, beta, &
1821 : C_re, C_im, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, &
1822 : c_first_row)
1823 : CHARACTER(LEN=1), INTENT(IN) :: transa, transb
1824 : INTEGER, INTENT(IN) :: m, n, k
1825 : REAL(KIND=dp), INTENT(IN) :: alpha
1826 : TYPE(cp_fm_type), INTENT(IN) :: A_re, A_im, B_re, B_im
1827 : REAL(KIND=dp), INTENT(IN) :: beta
1828 : TYPE(cp_fm_type), INTENT(IN) :: C_re, C_im
1829 : INTEGER, INTENT(IN), OPTIONAL :: a_first_col, a_first_row, b_first_col, &
1830 : b_first_row, c_first_col, c_first_row
1831 :
1832 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_complex_fm_gemm'
1833 :
1834 : CHARACTER(LEN=1) :: transa_real, transb_real
1835 : INTEGER :: handle
1836 : REAL(KIND=dp) :: sign_a, sign_b
1837 :
1838 4748 : CALL timeset(routineN, handle)
1839 :
1840 4748 : transa_real = transa
1841 4748 : transb_real = transb
1842 4748 : sign_a = 1.0_dp
1843 4748 : sign_b = 1.0_dp
1844 4748 : IF (transa == "C") THEN
1845 3364 : transa_real = "T"
1846 3364 : sign_a = -1.0_dp
1847 : END IF
1848 4748 : IF (transb == "C") THEN
1849 240 : transb_real = "T"
1850 240 : sign_b = -1.0_dp
1851 : END IF
1852 :
1853 : CALL cp_fm_gemm(transa_real, transb_real, m, n, k, alpha, A_re, B_re, beta, C_re, &
1854 : a_first_col=a_first_col, &
1855 : a_first_row=a_first_row, &
1856 : b_first_col=b_first_col, &
1857 : b_first_row=b_first_row, &
1858 : c_first_col=c_first_col, &
1859 4748 : c_first_row=c_first_row)
1860 : CALL cp_fm_gemm(transa_real, transb_real, m, n, k, -sign_a*sign_b*alpha, &
1861 : A_im, B_im, 1.0_dp, C_re, &
1862 : a_first_col=a_first_col, &
1863 : a_first_row=a_first_row, &
1864 : b_first_col=b_first_col, &
1865 : b_first_row=b_first_row, &
1866 : c_first_col=c_first_col, &
1867 4748 : c_first_row=c_first_row)
1868 : CALL cp_fm_gemm(transa_real, transb_real, m, n, k, sign_b*alpha, A_re, B_im, beta, C_im, &
1869 : a_first_col=a_first_col, &
1870 : a_first_row=a_first_row, &
1871 : b_first_col=b_first_col, &
1872 : b_first_row=b_first_row, &
1873 : c_first_col=c_first_col, &
1874 4748 : c_first_row=c_first_row)
1875 : CALL cp_fm_gemm(transa_real, transb_real, m, n, k, sign_a*alpha, A_im, B_re, 1.0_dp, C_im, &
1876 : a_first_col=a_first_col, &
1877 : a_first_row=a_first_row, &
1878 : b_first_col=b_first_col, &
1879 : b_first_row=b_first_row, &
1880 : c_first_col=c_first_col, &
1881 4748 : c_first_row=c_first_row)
1882 :
1883 4748 : CALL timestop(handle)
1884 :
1885 4748 : END SUBROUTINE cp_complex_fm_gemm
1886 :
1887 : ! **************************************************************************************************
1888 : !> \brief inverts a matrix using LU decomposition
1889 : !> the input matrix will be overwritten
1890 : !> \param matrix : input a general square non-singular matrix, outputs its inverse
1891 : !> \param info_out : optional, if present outputs the info from (p)zgetri
1892 : !> \author Lianheng Tong
1893 : ! **************************************************************************************************
1894 0 : SUBROUTINE cp_fm_lu_invert(matrix, info_out)
1895 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1896 : INTEGER, INTENT(OUT), OPTIONAL :: info_out
1897 :
1898 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_lu_invert'
1899 :
1900 : INTEGER :: nrows_global, handle, info, lwork
1901 0 : INTEGER, DIMENSION(:), ALLOCATABLE :: ipivot
1902 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mat
1903 0 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
1904 : #if defined(__parallel)
1905 : INTEGER :: liwork
1906 : INTEGER, DIMENSION(9) :: desca
1907 0 : INTEGER, DIMENSION(:), ALLOCATABLE :: iwork
1908 : #else
1909 : INTEGER :: lda
1910 : #endif
1911 :
1912 0 : CALL timeset(routineN, handle)
1913 :
1914 0 : mat => matrix%local_data
1915 0 : nrows_global = matrix%matrix_struct%nrow_global
1916 0 : CPASSERT(nrows_global == matrix%matrix_struct%ncol_global)
1917 0 : ALLOCATE (ipivot(nrows_global))
1918 : ! do LU decomposition
1919 : #if defined(__parallel)
1920 0 : desca = matrix%matrix_struct%descriptor
1921 : CALL pdgetrf(nrows_global, nrows_global, &
1922 0 : mat, 1, 1, desca, ipivot, info)
1923 : #else
1924 : lda = SIZE(mat, 1)
1925 : CALL dgetrf(nrows_global, nrows_global, &
1926 : mat, lda, ipivot, info)
1927 : #endif
1928 0 : IF (info /= 0) THEN
1929 0 : CALL cp_abort(__LOCATION__, "LU decomposition has failed")
1930 : END IF
1931 : ! do inversion
1932 : #if defined(__parallel)
1933 0 : ALLOCATE (iwork(1))
1934 : CALL pdgetri(nrows_global, mat, 1, 1, desca, &
1935 0 : ipivot, work, -1, iwork, -1, info)
1936 : lwork = INT(work(1))
1937 0 : DEALLOCATE (work)
1938 : ALLOCATE (work(lwork))
1939 : liwork = INT(iwork(1))
1940 : DEALLOCATE (iwork)
1941 : ALLOCATE (iwork(liwork))
1942 : CALL pdgetri(nrows_global, mat, 1, 1, desca, &
1943 : ipivot, work, lwork, iwork, liwork, info)
1944 : DEALLOCATE (iwork)
1945 : #else
1946 : CALL dgetri(nrows_global, mat, lda, &
1947 : ipivot, work, -1, info)
1948 : lwork = INT(work(1))
1949 : DEALLOCATE (work)
1950 : ALLOCATE (work(lwork))
1951 : CALL dgetri(nrows_global, mat, lda, &
1952 : ipivot, work, lwork, info)
1953 : #endif
1954 : DEALLOCATE (work)
1955 : DEALLOCATE (ipivot)
1956 :
1957 : IF (PRESENT(info_out)) THEN
1958 : info_out = info
1959 : ELSE
1960 : IF (info /= 0) &
1961 : CALL cp_abort(__LOCATION__, "LU inversion has failed")
1962 : END IF
1963 :
1964 : CALL timestop(handle)
1965 :
1966 0 : END SUBROUTINE cp_fm_lu_invert
1967 :
1968 : ! **************************************************************************************************
1969 : !> \brief norm of matrix using (p)dlange
1970 : !> \param matrix : input a general matrix
1971 : !> \param mode : 'M' max abs element value,
1972 : !> '1' or 'O' one norm, i.e. maximum column sum
1973 : !> 'I' infinity norm, i.e. maximum row sum
1974 : !> 'F' or 'E' Frobenius norm, i.e. sqrt of sum of all squares of elements
1975 : !> \return : the norm according to mode
1976 : !> \author Lianheng Tong
1977 : ! **************************************************************************************************
1978 1332 : FUNCTION cp_fm_norm(matrix, mode) RESULT(res)
1979 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1980 : CHARACTER, INTENT(IN) :: mode
1981 : REAL(KIND=dp) :: res
1982 :
1983 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_norm'
1984 :
1985 : INTEGER :: nrows, ncols, handle, lwork, nrows_local, ncols_local
1986 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: aa
1987 1332 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
1988 : #if defined(__parallel)
1989 : INTEGER, DIMENSION(9) :: desca
1990 : #else
1991 : INTEGER :: lda
1992 : #endif
1993 :
1994 1332 : CALL timeset(routineN, handle)
1995 :
1996 : CALL cp_fm_get_info(matrix=matrix, &
1997 : nrow_global=nrows, &
1998 : ncol_global=ncols, &
1999 : nrow_local=nrows_local, &
2000 1332 : ncol_local=ncols_local)
2001 1332 : aa => matrix%local_data
2002 :
2003 : #if defined(__parallel)
2004 13320 : desca = matrix%matrix_struct%descriptor
2005 : SELECT CASE (mode)
2006 : CASE ('M', 'm')
2007 1332 : lwork = 1
2008 : CASE ('1', 'O', 'o')
2009 1332 : lwork = ncols_local
2010 : CASE ('I', 'i')
2011 0 : lwork = nrows_local
2012 : CASE ('F', 'f', 'E', 'e')
2013 0 : lwork = 1
2014 : CASE DEFAULT
2015 1332 : CPABORT("mode input is not valid")
2016 : END SELECT
2017 3996 : ALLOCATE (work(lwork))
2018 1332 : res = pdlange(mode, nrows, ncols, aa, 1, 1, desca, work)
2019 1332 : DEALLOCATE (work)
2020 : #else
2021 : SELECT CASE (mode)
2022 : CASE ('M', 'm')
2023 : lwork = 1
2024 : CASE ('1', 'O', 'o')
2025 : lwork = 1
2026 : CASE ('I', 'i')
2027 : lwork = nrows
2028 : CASE ('F', 'f', 'E', 'e')
2029 : lwork = 1
2030 : CASE DEFAULT
2031 : CPABORT("mode input is not valid")
2032 : END SELECT
2033 : ALLOCATE (work(lwork))
2034 : lda = SIZE(aa, 1)
2035 : res = dlange(mode, nrows, ncols, aa, lda, work)
2036 : DEALLOCATE (work)
2037 : #endif
2038 :
2039 1332 : CALL timestop(handle)
2040 :
2041 1332 : END FUNCTION cp_fm_norm
2042 :
2043 : ! **************************************************************************************************
2044 : !> \brief trace of a matrix using pdlatra
2045 : !> \param matrix : input a square matrix
2046 : !> \return : the trace
2047 : !> \author Lianheng Tong
2048 : ! **************************************************************************************************
2049 0 : FUNCTION cp_fm_latra(matrix) RESULT(res)
2050 : TYPE(cp_fm_type), INTENT(IN) :: matrix
2051 : REAL(KIND=dp) :: res
2052 :
2053 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_latra'
2054 :
2055 : INTEGER :: nrows, ncols, handle
2056 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: aa
2057 : #if defined(__parallel)
2058 : INTEGER, DIMENSION(9) :: desca
2059 : #else
2060 : INTEGER :: ii
2061 : #endif
2062 :
2063 0 : CALL timeset(routineN, handle)
2064 :
2065 0 : nrows = matrix%matrix_struct%nrow_global
2066 0 : ncols = matrix%matrix_struct%ncol_global
2067 0 : CPASSERT(nrows == ncols)
2068 0 : aa => matrix%local_data
2069 :
2070 : #if defined(__parallel)
2071 0 : desca = matrix%matrix_struct%descriptor
2072 0 : res = pdlatra(nrows, aa, 1, 1, desca)
2073 : #else
2074 : res = 0.0_dp
2075 : DO ii = 1, nrows
2076 : res = res + aa(ii, ii)
2077 : END DO
2078 : #endif
2079 :
2080 0 : CALL timestop(handle)
2081 :
2082 0 : END FUNCTION cp_fm_latra
2083 :
2084 : ! **************************************************************************************************
2085 : !> \brief compute a QR factorization with column pivoting of a M-by-N distributed matrix
2086 : !> sub( A ) = A(IA:IA+M-1,JA:JA+N-1)
2087 : !> \param matrix : input M-by-N distributed matrix sub( A ) which is to be factored
2088 : !> \param tau : scalar factors TAU of the elementary reflectors. TAU is tied to the distributed matrix A
2089 : !> \param nrow ...
2090 : !> \param ncol ...
2091 : !> \param first_row ...
2092 : !> \param first_col ...
2093 : !> \author MI
2094 : ! **************************************************************************************************
2095 38 : SUBROUTINE cp_fm_pdgeqpf(matrix, tau, nrow, ncol, first_row, first_col)
2096 :
2097 : TYPE(cp_fm_type), INTENT(IN) :: matrix
2098 : REAL(KIND=dp), DIMENSION(:), POINTER :: tau
2099 : INTEGER, INTENT(IN) :: nrow, ncol, first_row, first_col
2100 :
2101 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_pdgeqpf'
2102 :
2103 : INTEGER :: handle
2104 : INTEGER :: info, lwork
2105 38 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
2106 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
2107 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
2108 : #if defined(__parallel)
2109 : INTEGER, DIMENSION(9) :: descc
2110 : #else
2111 : INTEGER :: lda
2112 : #endif
2113 :
2114 38 : CALL timeset(routineN, handle)
2115 :
2116 38 : a => matrix%local_data
2117 38 : lwork = -1
2118 114 : ALLOCATE (work(2*nrow))
2119 114 : ALLOCATE (ipiv(ncol))
2120 38 : info = 0
2121 :
2122 : #if defined(__parallel)
2123 380 : descc(:) = matrix%matrix_struct%descriptor(:)
2124 : ! Call SCALAPACK routine to get optimal work dimension
2125 38 : CALL pdgeqpf(nrow, ncol, a, first_row, first_col, descc, ipiv, tau, work, lwork, info)
2126 38 : lwork = INT(work(1))
2127 38 : DEALLOCATE (work)
2128 114 : ALLOCATE (work(lwork))
2129 254 : tau = 0.0_dp
2130 38 : ipiv = 0
2131 :
2132 : ! Call SCALAPACK routine to get QR decomposition of CTs
2133 38 : CALL pdgeqpf(nrow, ncol, a, first_row, first_col, descc, ipiv, tau, work, lwork, info)
2134 : #else
2135 : CPASSERT(first_row == 1 .AND. first_col == 1)
2136 : lda = SIZE(a, 1)
2137 : CALL dgeqp3(nrow, ncol, a, lda, ipiv, tau, work, lwork, info)
2138 : lwork = INT(work(1))
2139 : DEALLOCATE (work)
2140 : ALLOCATE (work(lwork))
2141 : tau = 0.0_dp
2142 : ipiv = 0
2143 : CALL dgeqp3(nrow, ncol, a, lda, ipiv, tau, work, lwork, info)
2144 : #endif
2145 38 : CPASSERT(info == 0)
2146 :
2147 38 : DEALLOCATE (work)
2148 38 : DEALLOCATE (ipiv)
2149 :
2150 38 : CALL timestop(handle)
2151 :
2152 38 : END SUBROUTINE cp_fm_pdgeqpf
2153 :
2154 : ! **************************************************************************************************
2155 : !> \brief generates an M-by-N real distributed matrix Q denoting A(IA:IA+M-1,JA:JA+N-1)
2156 : !> with orthonormal columns, which is defined as the first N columns of a product of K
2157 : !> elementary reflectors of order M
2158 : !> \param matrix : On entry, the j-th column must contain the vector which defines the elementary reflector
2159 : !> as returned from PDGEQRF
2160 : !> On exit it contains the M-by-N distributed matrix Q
2161 : !> \param tau : contains the scalar factors TAU of elementary reflectors as returned by PDGEQRF
2162 : !> \param nrow ...
2163 : !> \param first_row ...
2164 : !> \param first_col ...
2165 : !> \author MI
2166 : ! **************************************************************************************************
2167 38 : SUBROUTINE cp_fm_pdorgqr(matrix, tau, nrow, first_row, first_col)
2168 :
2169 : TYPE(cp_fm_type), INTENT(IN) :: matrix
2170 : REAL(KIND=dp), DIMENSION(:), POINTER :: tau
2171 : INTEGER, INTENT(IN) :: nrow, first_row, first_col
2172 :
2173 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_pdorgqr'
2174 :
2175 : INTEGER :: handle
2176 : INTEGER :: info, lwork
2177 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
2178 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
2179 : #if defined(__parallel)
2180 : INTEGER, DIMENSION(9) :: descc
2181 : #else
2182 : INTEGER :: lda
2183 : #endif
2184 :
2185 38 : CALL timeset(routineN, handle)
2186 :
2187 38 : a => matrix%local_data
2188 38 : lwork = -1
2189 114 : ALLOCATE (work(2*nrow))
2190 38 : info = 0
2191 :
2192 : #if defined(__parallel)
2193 380 : descc(:) = matrix%matrix_struct%descriptor(:)
2194 :
2195 38 : CALL pdorgqr(nrow, nrow, nrow, a, first_row, first_col, descc, tau, work, lwork, info)
2196 38 : CPASSERT(info == 0)
2197 38 : lwork = INT(work(1))
2198 38 : DEALLOCATE (work)
2199 114 : ALLOCATE (work(lwork))
2200 :
2201 : ! Call SCALAPACK routine to get Q
2202 38 : CALL pdorgqr(nrow, nrow, nrow, a, first_row, first_col, descc, tau, work, lwork, info)
2203 : #else
2204 : CPASSERT(first_row == 1 .AND. first_col == 1)
2205 : lda = SIZE(a, 1)
2206 : CALL dorgqr(nrow, nrow, nrow, a, lda, tau, work, lwork, info)
2207 : lwork = INT(work(1))
2208 : DEALLOCATE (work)
2209 : ALLOCATE (work(lwork))
2210 : CALL dorgqr(nrow, nrow, nrow, a, lda, tau, work, lwork, info)
2211 : #endif
2212 38 : CPASSERT(INFO == 0)
2213 :
2214 38 : DEALLOCATE (work)
2215 38 : CALL timestop(handle)
2216 :
2217 38 : END SUBROUTINE cp_fm_pdorgqr
2218 :
2219 : ! **************************************************************************************************
2220 : !> \brief Applies a planar rotation defined by cs and sn to the i'th and j'th rows.
2221 : !> \param matrix ...
2222 : !> \param irow ...
2223 : !> \param jrow ...
2224 : !> \param cs cosine of the rotation angle
2225 : !> \param sn sinus of the rotation angle
2226 : !> \author Ole Schuett
2227 : ! **************************************************************************************************
2228 31928 : SUBROUTINE cp_fm_rot_rows(matrix, irow, jrow, cs, sn)
2229 : TYPE(cp_fm_type), INTENT(IN) :: matrix
2230 : INTEGER, INTENT(IN) :: irow, jrow
2231 : REAL(dp), INTENT(IN) :: cs, sn
2232 :
2233 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_rot_rows'
2234 : INTEGER :: handle, ncol
2235 :
2236 : #if defined(__parallel)
2237 : INTEGER :: info, lwork
2238 : INTEGER, DIMENSION(9) :: desc
2239 31928 : REAL(dp), DIMENSION(:), ALLOCATABLE :: work
2240 : #endif
2241 31928 : CALL timeset(routineN, handle)
2242 31928 : CALL cp_fm_get_info(matrix, ncol_global=ncol)
2243 : #if defined(__parallel)
2244 31928 : IF (1 /= matrix%matrix_struct%context%n_pid) THEN
2245 31928 : lwork = 2*ncol + 1
2246 95784 : ALLOCATE (work(lwork))
2247 319280 : desc(:) = matrix%matrix_struct%descriptor(:)
2248 : CALL pdrot(ncol, &
2249 : matrix%local_data(1, 1), irow, 1, desc, ncol, &
2250 : matrix%local_data(1, 1), jrow, 1, desc, ncol, &
2251 31928 : cs, sn, work, lwork, info)
2252 31928 : CPASSERT(info == 0)
2253 31928 : DEALLOCATE (work)
2254 : ELSE
2255 : #endif
2256 0 : CALL drot(ncol, matrix%local_data(irow, 1), ncol, matrix%local_data(jrow, 1), ncol, cs, sn)
2257 : #if defined(__parallel)
2258 : END IF
2259 : #endif
2260 31928 : CALL timestop(handle)
2261 31928 : END SUBROUTINE cp_fm_rot_rows
2262 :
2263 : ! **************************************************************************************************
2264 : !> \brief Applies a planar rotation defined by cs and sn to the i'th and j'th columnns.
2265 : !> \param matrix ...
2266 : !> \param icol ...
2267 : !> \param jcol ...
2268 : !> \param cs cosine of the rotation angle
2269 : !> \param sn sinus of the rotation angle
2270 : !> \author Ole Schuett
2271 : ! **************************************************************************************************
2272 36928 : SUBROUTINE cp_fm_rot_cols(matrix, icol, jcol, cs, sn)
2273 : TYPE(cp_fm_type), INTENT(IN) :: matrix
2274 : INTEGER, INTENT(IN) :: icol, jcol
2275 : REAL(dp), INTENT(IN) :: cs, sn
2276 :
2277 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_rot_cols'
2278 : INTEGER :: handle, nrow
2279 :
2280 : #if defined(__parallel)
2281 : INTEGER :: info, lwork
2282 : INTEGER, DIMENSION(9) :: desc
2283 36928 : REAL(dp), DIMENSION(:), ALLOCATABLE :: work
2284 : #endif
2285 36928 : CALL timeset(routineN, handle)
2286 36928 : CALL cp_fm_get_info(matrix, nrow_global=nrow)
2287 : #if defined(__parallel)
2288 36928 : IF (1 /= matrix%matrix_struct%context%n_pid) THEN
2289 36928 : lwork = 2*nrow + 1
2290 110784 : ALLOCATE (work(lwork))
2291 369280 : desc(:) = matrix%matrix_struct%descriptor(:)
2292 : CALL pdrot(nrow, &
2293 : matrix%local_data(1, 1), 1, icol, desc, 1, &
2294 : matrix%local_data(1, 1), 1, jcol, desc, 1, &
2295 36928 : cs, sn, work, lwork, info)
2296 36928 : CPASSERT(info == 0)
2297 36928 : DEALLOCATE (work)
2298 : ELSE
2299 : #endif
2300 0 : CALL drot(nrow, matrix%local_data(1, icol), 1, matrix%local_data(1, jcol), 1, cs, sn)
2301 : #if defined(__parallel)
2302 : END IF
2303 : #endif
2304 36928 : CALL timestop(handle)
2305 36928 : END SUBROUTINE cp_fm_rot_cols
2306 :
2307 : ! **************************************************************************************************
2308 : !> \brief Orthonormalizes selected rows and columns of a full matrix, matrix_a
2309 : !> \param matrix_a ...
2310 : !> \param B ...
2311 : !> \param nrows number of rows of matrix_a, optional, defaults to size(matrix_a,1)
2312 : !> \param ncols number of columns of matrix_a, optional, defaults to size(matrix_a, 2)
2313 : !> \param start_row starting index of rows, optional, defaults to 1
2314 : !> \param start_col starting index of columns, optional, defaults to 1
2315 : !> \param do_norm ...
2316 : !> \param do_print ...
2317 : ! **************************************************************************************************
2318 0 : SUBROUTINE cp_fm_Gram_Schmidt_orthonorm(matrix_a, B, nrows, ncols, start_row, start_col, &
2319 : do_norm, do_print)
2320 :
2321 : TYPE(cp_fm_type), INTENT(IN) :: matrix_a
2322 : REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: B
2323 : INTEGER, INTENT(IN), OPTIONAL :: nrows, ncols, start_row, start_col
2324 : LOGICAL, INTENT(IN), OPTIONAL :: do_norm, do_print
2325 :
2326 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_Gram_Schmidt_orthonorm'
2327 :
2328 : INTEGER :: end_col_global, end_col_local, end_row_global, end_row_local, handle, i, j, &
2329 : j_col, ncol_global, ncol_local, nrow_global, nrow_local, start_col_global, &
2330 : start_col_local, start_row_global, start_row_local, this_col, unit_nr
2331 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2332 : LOGICAL :: my_do_norm, my_do_print
2333 : REAL(KIND=dp) :: norm
2334 0 : REAL(kind=dp), DIMENSION(:, :), POINTER :: a
2335 :
2336 0 : CALL timeset(routineN, handle)
2337 :
2338 0 : my_do_norm = .TRUE.
2339 0 : IF (PRESENT(do_norm)) my_do_norm = do_norm
2340 :
2341 0 : my_do_print = .FALSE.
2342 0 : IF (PRESENT(do_print) .AND. (my_do_norm)) my_do_print = do_print
2343 :
2344 0 : unit_nr = -1
2345 0 : IF (my_do_print) THEN
2346 0 : unit_nr = cp_logger_get_default_unit_nr()
2347 0 : IF (unit_nr < 1) my_do_print = .FALSE.
2348 : END IF
2349 :
2350 0 : IF (SIZE(B) /= 0) THEN
2351 0 : IF (PRESENT(nrows)) THEN
2352 0 : nrow_global = nrows
2353 : ELSE
2354 0 : nrow_global = SIZE(B, 1)
2355 : END IF
2356 :
2357 0 : IF (PRESENT(ncols)) THEN
2358 0 : ncol_global = ncols
2359 : ELSE
2360 0 : ncol_global = SIZE(B, 2)
2361 : END IF
2362 :
2363 0 : IF (PRESENT(start_row)) THEN
2364 0 : start_row_global = start_row
2365 : ELSE
2366 : start_row_global = 1
2367 : END IF
2368 :
2369 0 : IF (PRESENT(start_col)) THEN
2370 0 : start_col_global = start_col
2371 : ELSE
2372 : start_col_global = 1
2373 : END IF
2374 :
2375 0 : end_row_global = start_row_global + nrow_global - 1
2376 0 : end_col_global = start_col_global + ncol_global - 1
2377 :
2378 : CALL cp_fm_get_info(matrix=matrix_a, &
2379 : nrow_global=nrow_global, ncol_global=ncol_global, &
2380 : nrow_local=nrow_local, ncol_local=ncol_local, &
2381 0 : row_indices=row_indices, col_indices=col_indices)
2382 0 : IF (end_row_global > nrow_global) THEN
2383 : end_row_global = nrow_global
2384 : END IF
2385 0 : IF (end_col_global > ncol_global) THEN
2386 : end_col_global = ncol_global
2387 : END IF
2388 :
2389 : ! find out row/column indices of locally stored matrix elements that
2390 : ! needs to be copied.
2391 : ! Arrays row_indices and col_indices are assumed to be sorted in
2392 : ! ascending order
2393 0 : DO start_row_local = 1, nrow_local
2394 0 : IF (row_indices(start_row_local) >= start_row_global) EXIT
2395 : END DO
2396 :
2397 0 : DO end_row_local = start_row_local, nrow_local
2398 0 : IF (row_indices(end_row_local) > end_row_global) EXIT
2399 : END DO
2400 0 : end_row_local = end_row_local - 1
2401 :
2402 0 : DO start_col_local = 1, ncol_local
2403 0 : IF (col_indices(start_col_local) >= start_col_global) EXIT
2404 : END DO
2405 :
2406 0 : DO end_col_local = start_col_local, ncol_local
2407 0 : IF (col_indices(end_col_local) > end_col_global) EXIT
2408 : END DO
2409 0 : end_col_local = end_col_local - 1
2410 :
2411 0 : a => matrix_a%local_data
2412 :
2413 0 : this_col = col_indices(start_col_local) - start_col_global + 1
2414 :
2415 0 : B(:, this_col) = a(:, start_col_local)
2416 :
2417 0 : IF (my_do_norm) THEN
2418 0 : norm = SQRT(accurate_dot_product(B(:, this_col), B(:, this_col)))
2419 0 : B(:, this_col) = B(:, this_col)/norm
2420 0 : IF (my_do_print) WRITE (unit_nr, '(I3,F8.3)') this_col, norm
2421 : END IF
2422 :
2423 0 : DO i = start_col_local + 1, end_col_local
2424 0 : this_col = col_indices(i) - start_col_global + 1
2425 0 : B(:, this_col) = a(:, i)
2426 0 : DO j = start_col_local, i - 1
2427 0 : j_col = col_indices(j) - start_col_global + 1
2428 : B(:, this_col) = B(:, this_col) - &
2429 : accurate_dot_product(B(:, j_col), B(:, this_col))* &
2430 0 : B(:, j_col)/accurate_dot_product(B(:, j_col), B(:, j_col))
2431 : END DO
2432 :
2433 0 : IF (my_do_norm) THEN
2434 0 : norm = SQRT(accurate_dot_product(B(:, this_col), B(:, this_col)))
2435 0 : B(:, this_col) = B(:, this_col)/norm
2436 0 : IF (my_do_print) WRITE (unit_nr, '(I3,F8.3)') this_col, norm
2437 : END IF
2438 :
2439 : END DO
2440 0 : CALL matrix_a%matrix_struct%para_env%sum(B)
2441 : END IF
2442 :
2443 0 : CALL timestop(handle)
2444 :
2445 0 : END SUBROUTINE cp_fm_Gram_Schmidt_orthonorm
2446 :
2447 : ! **************************************************************************************************
2448 : !> \brief Cholesky decomposition
2449 : !> \param fm_matrix ...
2450 : !> \param n ...
2451 : !> \param uplo triangular format; defaults to 'U'
2452 : ! **************************************************************************************************
2453 11255 : SUBROUTINE cp_fm_potrf(fm_matrix, n, uplo)
2454 : TYPE(cp_fm_type) :: fm_matrix
2455 : INTEGER, INTENT(IN) :: n
2456 : CHARACTER, INTENT(IN), OPTIONAL :: uplo
2457 :
2458 : CHARACTER :: myuplo
2459 : INTEGER :: info
2460 11255 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
2461 : #if defined(__parallel)
2462 : INTEGER, DIMENSION(9) :: desca
2463 : #endif
2464 :
2465 11255 : myuplo = 'U'
2466 0 : IF (PRESENT(uplo)) myuplo = uplo
2467 :
2468 11255 : a => fm_matrix%local_data
2469 : #if defined(__parallel)
2470 112550 : desca(:) = fm_matrix%matrix_struct%descriptor(:)
2471 11255 : CALL pdpotrf(myuplo, n, a(1, 1), 1, 1, desca, info)
2472 : #else
2473 : CALL dpotrf(myuplo, n, a(1, 1), SIZE(a, 1), info)
2474 : #endif
2475 11255 : IF (info /= 0) &
2476 0 : CPABORT("Cholesky decomposition failed. Matrix ill-conditioned?")
2477 :
2478 11255 : END SUBROUTINE cp_fm_potrf
2479 :
2480 : ! **************************************************************************************************
2481 : !> \brief Invert trianguar matrix
2482 : !> \param fm_matrix the matrix to invert (triangular matrix according to uplo)
2483 : !> \param n size of the matrix to invert
2484 : !> \param uplo triangular format; defaults to 'U'
2485 : ! **************************************************************************************************
2486 10432 : SUBROUTINE cp_fm_potri(fm_matrix, n, uplo)
2487 : TYPE(cp_fm_type) :: fm_matrix
2488 : INTEGER, INTENT(IN) :: n
2489 : CHARACTER, INTENT(IN), OPTIONAL :: uplo
2490 :
2491 : CHARACTER :: myuplo
2492 10432 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
2493 : INTEGER :: info
2494 : #if defined(__parallel)
2495 : INTEGER, DIMENSION(9) :: desca
2496 : #endif
2497 :
2498 10432 : myuplo = 'U'
2499 0 : IF (PRESENT(uplo)) myuplo = uplo
2500 :
2501 10432 : a => fm_matrix%local_data
2502 : #if defined(__parallel)
2503 104320 : desca(:) = fm_matrix%matrix_struct%descriptor(:)
2504 10432 : CALL pdpotri(myuplo, n, a(1, 1), 1, 1, desca, info)
2505 : #else
2506 : CALL dpotri(myuplo, n, a(1, 1), SIZE(a, 1), info)
2507 : #endif
2508 10432 : CPASSERT(info == 0)
2509 10432 : END SUBROUTINE cp_fm_potri
2510 :
2511 : ! **************************************************************************************************
2512 : !> \brief Calculates
2513 : !> yv = alpha*amat*xv + beta*yv
2514 : !> where amat: fm matrix
2515 : !> xv : vector replicated
2516 : !> yv : vector replicated
2517 : !> Defaults: alpha = 1, beta = 0
2518 : ! **************************************************************************************************
2519 22982 : SUBROUTINE cp_fm_matvec(amat, xv, yv, alpha, beta)
2520 : TYPE(cp_fm_type), INTENT(IN) :: amat
2521 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: xv
2522 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: yv
2523 : REAL(KIND=dp), OPTIONAL, INTENT(IN) :: alpha, beta
2524 :
2525 : INTEGER :: na, nc, nx, ny
2526 : REAL(KIND=dp) :: aval, bval
2527 : #if defined(__parallel)
2528 : INTEGER :: nrl, ncl, ic, ir
2529 22982 : INTEGER, DIMENSION(:), POINTER :: rind, cind
2530 22982 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: xvl, yvl, yvm
2531 : #endif
2532 :
2533 22982 : aval = 1.0_dp
2534 22480 : IF (PRESENT(alpha)) aval = alpha
2535 22982 : bval = 0.0_dp
2536 22982 : IF (PRESENT(beta)) bval = beta
2537 :
2538 22982 : CALL cp_fm_get_info(amat, nrow_global=na, ncol_global=nc)
2539 22982 : nx = SIZE(xv)
2540 22982 : ny = SIZE(yv)
2541 22982 : IF ((nx /= ny) .OR. (nc /= nx)) THEN
2542 0 : CPABORT("cp_fm_matvec: incompatible dimensions")
2543 : END IF
2544 : #if defined(__parallel)
2545 : CALL cp_fm_get_info(amat, nrow_local=nrl, ncol_local=ncl, &
2546 22982 : row_indices=rind, col_indices=cind)
2547 160874 : ALLOCATE (xvl(ncl), yvl(nrl), yvm(ny))
2548 375484 : DO ic = 1, ncl
2549 375484 : xvl(ic) = xv(cind(ic))
2550 : END DO
2551 22982 : yvl = 0.0_dp
2552 22982 : IF (nrl > 0 .AND. ncl > 0) THEN
2553 15462576 : yvl(1:nrl) = MATMUL(amat%local_data(1:nrl, 1:ncl), xvl)
2554 : END IF
2555 22982 : yvm = 0.0_dp
2556 199233 : DO ir = 1, nrl
2557 199233 : yvm(rind(ir)) = yvl(ir)
2558 : END DO
2559 22982 : CALL amat%matrix_struct%para_env%sum(yvm)
2560 22982 : IF (bval == 0.0_dp) THEN
2561 280306 : yv = aval*yvm
2562 : ELSE
2563 95178 : yv = bval*yv + aval*yvm
2564 : END IF
2565 : #else
2566 : IF (bval == 0.0_dp) THEN
2567 : yv = aval*MATMUL(amat%local_data, xv)
2568 : ELSE
2569 : yv = bval*yv + aval*MATMUL(amat%local_data, xv)
2570 : END IF
2571 : #endif
2572 :
2573 91928 : END SUBROUTINE cp_fm_matvec
2574 :
2575 : ! **************************************************************************************************
2576 : !> \brief Adds a scalar to the diagonal of a distributed real full matrix.
2577 : !> \param matrix input/output matrix
2578 : !> \param alpha scalar added to each diagonal element
2579 : !> \param n_active optional number of physical rows/columns
2580 : ! **************************************************************************************************
2581 53946 : SUBROUTINE cp_fm_add_on_diag(matrix, alpha, n_active)
2582 : TYPE(cp_fm_type), INTENT(IN) :: matrix
2583 : REAL(KIND=dp), INTENT(IN) :: alpha
2584 : INTEGER, INTENT(IN), OPTIONAL :: n_active
2585 :
2586 : INTEGER :: i_global, n_diag, ncol_global, nrow_global
2587 : #if defined(__parallel)
2588 : INTEGER :: icol_local, ipcol, iprow, irow, irow_local, &
2589 : mypcol, myprow, npcol, nprow, nrow_local
2590 : INTEGER, DIMENSION(9) :: desca
2591 26973 : INTEGER, DIMENSION(:), POINTER :: row_indices
2592 : TYPE(cp_blacs_env_type), POINTER :: context
2593 : #endif
2594 :
2595 26973 : CALL cp_fm_get_info(matrix, nrow_global=nrow_global, ncol_global=ncol_global)
2596 :
2597 26973 : n_diag = MIN(nrow_global, ncol_global)
2598 26973 : IF (PRESENT(n_active)) n_diag = n_active
2599 26973 : CPASSERT(n_diag >= 0 .AND. n_diag <= MIN(nrow_global, ncol_global))
2600 :
2601 : #if defined(__parallel)
2602 26973 : CALL cp_fm_get_info(matrix, nrow_local=nrow_local, row_indices=row_indices)
2603 26973 : context => matrix%matrix_struct%context
2604 26973 : myprow = context%mepos(1)
2605 26973 : mypcol = context%mepos(2)
2606 26973 : nprow = context%num_pe(1)
2607 26973 : npcol = context%num_pe(2)
2608 269730 : desca(:) = matrix%matrix_struct%descriptor(:)
2609 :
2610 1169808 : DO irow = 1, nrow_local
2611 1142835 : i_global = row_indices(irow)
2612 1142835 : IF (i_global > n_diag) EXIT
2613 : CALL infog2l(i_global, i_global, desca, nprow, npcol, myprow, mypcol, &
2614 1142835 : irow_local, icol_local, iprow, ipcol)
2615 1169808 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
2616 1142227 : matrix%local_data(irow_local, icol_local) = matrix%local_data(irow_local, icol_local) + alpha
2617 : END IF
2618 : END DO
2619 : #else
2620 : DO i_global = 1, n_diag
2621 : matrix%local_data(i_global, i_global) = matrix%local_data(i_global, i_global) + alpha
2622 : END DO
2623 : #endif
2624 :
2625 26973 : END SUBROUTINE cp_fm_add_on_diag
2626 :
2627 : END MODULE cp_fm_basic_linalg
|
