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