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