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