1! 2! Copyright (C) 2016 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8MODULE matrix_inversion 9 10 IMPLICIT NONE 11 PRIVATE 12 PUBLIC :: invmat 13 14 INTERFACE invmat 15 MODULE PROCEDURE invmat_r, invmat_c 16 END INTERFACE 17 18 CONTAINS 19 20 SUBROUTINE invmat_r (n, a, a_inv, da) 21 !----------------------------------------------------------------------- 22 ! computes the inverse of a n*n real matrix "a" using LAPACK routines 23 ! if "a_inv" is not present, "a" contains the inverse on output 24 ! if "a_inv" is present, it contains the inverse on output, "a" is unchanged 25 ! if "da" is specified and if the matrix is dimensioned 3x3, 26 ! it also returns the determinant in "da" 27 ! 28 USE kinds, ONLY : DP 29 IMPLICIT NONE 30 INTEGER, INTENT(in) :: n 31 REAL(DP), DIMENSION (n,n), INTENT(inout) :: a 32 REAL(DP), DIMENSION (n,n), INTENT(out), OPTIONAL :: a_inv 33 REAL(DP), OPTIONAL, INTENT(out) :: da 34 ! 35 INTEGER :: info, lda, lwork 36 ! info=0: inversion was successful 37 ! lda : leading dimension (the same as n) 38 INTEGER, ALLOCATABLE :: ipiv (:) 39 ! ipiv : work space for pivoting 40 REAL(DP), ALLOCATABLE :: work (:) 41 ! more work space 42 INTEGER, SAVE :: lworkfact = 64 43 ! 44 IF ( PRESENT(da) ) THEN 45 IF ( n == 3 ) THEN 46 da = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) + & 47 a(1,2)*(a(2,3)*a(3,1)-a(2,1)*a(3,3)) + & 48 a(1,3)*(a(2,1)*a(3,2)-a(3,1)*a(2,2)) 49 IF (abs(da) < 1.d-10) CALL errore(' invmat ',' singular matrix ', 1) 50 ELSE 51 da = 0.0_dp 52 ENDIF 53 ENDIF 54 ! 55 lda = n 56 lwork=64*n 57 ALLOCATE(ipiv(n), work(lwork) ) 58 ! 59 IF ( PRESENT(a_inv) ) THEN 60 a_inv(:,:) = a(:,:) 61 CALL dgetrf (n, n, a_inv, lda, ipiv, info) 62 ELSE 63 CALL dgetrf (n, n, a, lda, ipiv, info) 64 END IF 65 CALL errore ('invmat', 'error in DGETRF', abs (info) ) 66 IF ( PRESENT(a_inv) ) THEN 67 CALL dgetri (n, a_inv, lda, ipiv, work, lwork, info) 68 ELSE 69 CALL dgetri (n, a, lda, ipiv, work, lwork, info) 70 END IF 71 CALL errore ('invmat', 'error in DGETRI', abs (info) ) 72 ! 73 lworkfact = INT (work(1)/n) 74 DEALLOCATE ( work, ipiv ) 75 76 END SUBROUTINE invmat_r 77 78 SUBROUTINE invmat_c (n, a, a_inv, da) 79 !----------------------------------------------------------------------- 80 ! as invmat_r, for a complex matrix 81 ! 82 USE kinds, ONLY : DP 83 IMPLICIT NONE 84 INTEGER, INTENT(IN) :: n 85 COMPLEX (DP), DIMENSION (n,n), INTENT(INOUT) :: a 86 COMPLEX (DP), OPTIONAL, DIMENSION (n,n), INTENT(OUT) :: a_inv 87 COMPLEX (DP), OPTIONAL, INTENT(OUT) :: da 88 ! 89 INTEGER :: info, lda, lwork 90 ! info=0: inversion was successful 91 ! lda : leading dimension (the same as n) 92 INTEGER, ALLOCATABLE :: ipiv (:) 93 ! ipiv : work space for pivoting 94 COMPLEX(DP), ALLOCATABLE :: work (:) 95 ! more work space 96 INTEGER, SAVE :: lworkfact = 64 97 ! 98 IF ( PRESENT(da) ) THEN 99 IF (n == 3) THEN 100 da = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) + & 101 a(1,2)*(a(2,3)*a(3,1)-a(2,1)*a(3,3)) + & 102 a(1,3)*(a(2,1)*a(3,2)-a(3,1)*a(2,2)) 103 IF (abs(da) < 1.d-10) CALL errore(' invmat ',' singular matrix ', 1) 104 ELSE 105 da = (0.d0,0.d0) 106 ENDIF 107 ENDIF 108 ! 109 lda = n 110 lwork=64*n 111 ALLOCATE(ipiv(n), work(lwork) ) 112 ! 113 IF ( PRESENT(a_inv) ) THEN 114 a_inv(:,:) = a(:,:) 115 CALL zgetrf (n, n, a_inv, lda, ipiv, info) 116 ELSE 117 CALL zgetrf (n, n, a, lda, ipiv, info) 118 END IF 119 CALL errore ('invmat', 'error in ZGETRF', abs (info) ) 120 IF ( PRESENT(a_inv) ) THEN 121 CALL zgetri (n, a_inv, lda, ipiv, work, lwork, info) 122 ELSE 123 CALL zgetri (n, a, lda, ipiv, work, lwork, info) 124 END IF 125 CALL errore ('invmat', 'error in ZGETRI', abs (info) ) 126 ! 127 lworkfact = INT (work(1)/n) 128 DEALLOCATE ( work, ipiv ) 129 ! 130 END SUBROUTINE invmat_c 131 132END MODULE matrix_inversion 133