types_rcmat.F90 6.62 KB
Newer Older
Daniel Wortmann's avatar
Daniel Wortmann committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14
module m_types_rcmat
  use m_judft
    IMPLICIT NONE
    TYPE t_mat
     LOGICAL :: l_real
     INTEGER :: matsize1=-1
     INTEGER :: matsize2=-1
     REAL,ALLOCATABLE    :: data_r(:,:)
     COMPLEX,ALLOCATABLE :: data_c(:,:)
   CONTAINS
     PROCEDURE        :: alloc => t_mat_alloc
     PROCEDURE        :: multiply=>t_mat_multiply
     PROCEDURE        :: transpose=>t_mat_transpose
     PROCEDURE        :: from_packed=>t_mat_from_packed
Daniel Wortmann's avatar
Daniel Wortmann committed
15 16
     PROCEDURE        :: inverse =>t_mat_inverse
     PROCEDURE        :: to_packed=>t_mat_to_packed
Daniel Wortmann's avatar
Daniel Wortmann committed
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
   END type t_mat

 CONTAINS
  SUBROUTINE t_mat_alloc(mat,l_real,matsize1,matsize2)
    CLASS(t_mat) :: mat
    LOGICAL,INTENT(IN),OPTIONAL:: l_real
    INTEGER,INTENT(IN),OPTIONAL:: matsize1,matsize2

    INTEGER:: err
    IF (present(l_real)) mat%l_real=l_real
    IF (present(matsize1)) mat%matsize1=matsize1
    IF (present(matsize2)) mat%matsize2=matsize2

    IF (mat%matsize1<0) CALL judft_error("Cannot allocate memory for mat datatype that is not initialized",hint="This is a BUG in FLEUR, please report")
    IF (mat%matsize2<0)  mat%matsize2=mat%matsize1 !Square by default
    
    IF (allocated(mat%data_r)) deallocate(mat%data_r)
    IF (allocated(mat%data_c)) deallocate(mat%data_c)
    
    IF (mat%l_real) THEN
       ALLOCATE(mat%data_r(mat%matsize1,mat%matsize2),STAT=err)
Daniel Wortmann's avatar
Daniel Wortmann committed
38
       ALLOCATE(mat%data_c(0,0))
Daniel Wortmann's avatar
Daniel Wortmann committed
39
    ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
40
       ALLOCATE(mat%data_r(0,0))
Daniel Wortmann's avatar
Daniel Wortmann committed
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
       ALLOCATE(mat%data_c(mat%matsize1,mat%matsize2),STAT=err)
    ENDIF

    IF (err>0) CALL judft_error("Allocation of memmory failed for mat datatype",hint="You probably run out of memory")
  END SUBROUTINE t_mat_alloc

  SUBROUTINE t_mat_multiply(mat1,mat2,res)
    CLASS(t_mat),INTENT(INOUT)       ::mat1
    TYPE(t_mat),INTENT(IN)           ::mat2
    TYPE(t_mat),INTENT(OUT),OPTIONAL ::res

    if (mat1%matsize2.ne.mat2%matsize1) CALL judft_error("Cannot multiply matrices because of non-matching dimensions",hint="This is a BUG in FLEUR, please report")
    
    IF (present(res)) THEN
       call res%alloc(mat1%l_real,mat1%matsize1,mat2%matsize2)
       IF (mat1%l_real) THEN
          res%data_r=matmul(mat1%data_r(:mat1%matsize1,:mat1%matsize2),mat2%data_r(:mat2%matsize1,:mat2%matsize2))
       ELSE
          res%data_c=matmul(mat1%data_c(:mat1%matsize1,:mat1%matsize2),mat2%data_c(:mat2%matsize1,:mat2%matsize2))
       ENDIF
    else
       if (mat1%matsize1.ne.mat1%matsize2) CALL judft_error("Cannot multiply matrices inplace because of non-matching dimensions",hint="This is a BUG in FLEUR, please report")
       if (mat1%l_real) THEN
          mat1%data_r(:mat1%matsize1,:mat1%matsize2)=matmul(mat1%data_r(:mat1%matsize1,:mat1%matsize2),mat2%data_r(:mat2%matsize1,:mat2%matsize2))
       ELSE
          mat1%data_c(:mat1%matsize1,:mat1%matsize2)=matmul(mat1%data_c(:mat1%matsize1,:mat1%matsize2),mat2%data_c(:mat2%matsize1,:mat2%matsize2))
       ENDIF
    end IF
  end SUBROUTINE t_mat_multiply


  SUBROUTINE t_mat_transpose(mat1,res)
    CLASS(t_mat),INTENT(INOUT)       ::mat1
    TYPE(t_mat),INTENT(OUT),OPTIONAL ::res
    
    IF (present(res)) THEN
       call res%alloc(mat1%l_real,mat1%matsize2,mat1%matsize1)
       IF (mat1%l_real) THEN
          res%data_r=transpose(mat1%data_r(:mat1%matsize1,:mat1%matsize2))
       ELSE
          res%data_c=transpose(mat1%data_c(:mat1%matsize1,:mat1%matsize2))
       ENDIF
    else
       if (mat1%matsize1.ne.mat1%matsize2) CALL judft_error("Cannot transpose matrices inplace because of non-matching dimensions",hint="This is a BUG in FLEUR, please report")
       IF (mat1%l_real) THEN
          mat1%data_r(:mat1%matsize1,:mat1%matsize2)=transpose(mat1%data_r(:mat1%matsize1,:mat1%matsize2))
       ELSE
          mat1%data_c(:mat1%matsize1,:mat1%matsize2)=transpose(mat1%data_c(:mat1%matsize1,:mat1%matsize2))
       ENDIF
    end IF
  end SUBROUTINE t_mat_transpose

Daniel Wortmann's avatar
Daniel Wortmann committed
93
  SUBROUTINE t_mat_from_packed(mat1,l_real,matsize,packed_r,packed_c)
Daniel Wortmann's avatar
Daniel Wortmann committed
94 95 96 97 98 99 100
    CLASS(t_mat),INTENT(INOUT)       :: mat1
    INTEGER,INTENT(IN)               :: matsize
    LOGICAL,INTENT(IN)               :: l_real
    REAL,INTENT(IN)                  :: packed_r(:)
    COMPLEX,INTENT(IN)               :: packed_c(:)

    INTEGER:: n,nn,i
Daniel Wortmann's avatar
Daniel Wortmann committed
101
    call mat1%alloc(l_real,matsize,matsize)
Daniel Wortmann's avatar
Daniel Wortmann committed
102 103 104 105 106
    i=1
    DO n=1,matsize
       DO nn=1,n
          if (l_real) THEN
             mat1%data_r(n,nn)=packed_r(i)
Daniel Wortmann's avatar
Daniel Wortmann committed
107
             mat1%data_r(nn,n)=packed_r(i)
Daniel Wortmann's avatar
Daniel Wortmann committed
108
          else
Daniel Wortmann's avatar
Daniel Wortmann committed
109 110
             mat1%data_c(n,nn)=conjg(packed_c(i))
             mat1%data_c(nn,n)=packed_c(i)
Daniel Wortmann's avatar
Daniel Wortmann committed
111 112 113 114 115
          end if
          i=i+1
       end DO
    end DO
  end SUBROUTINE t_mat_from_packed
Daniel Wortmann's avatar
Daniel Wortmann committed
116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164

  function t_mat_to_packed(mat)result(packed)
    CLASS(t_mat),INTENT(IN)       :: mat
    COMPLEX                       :: packed(mat%matsize1*(mat%matsize1+1)/2)
    integer :: n,nn,i
    real,parameter :: tol=1e-5
    if (mat%matsize1.ne.mat%matsize2) call judft_error("Could not pack no-square matrix",hint='This is a BUG, please report')
    i=1
    DO n=1,mat%matsize1
       DO nn=1,n
          if (mat%l_real) THEN
             packed(i)=(mat%data_r(n,nn)+mat%data_r(nn,n))/2.
             if (abs(mat%data_r(n,nn)-mat%data_r(nn,n))>tol) call judft_warn("Large unsymmetry in matrix packing")
          else
             packed(i)=(conjg(mat%data_c(n,nn))+mat%data_c(nn,n))/2.
             if (abs(conjg(mat%data_c(n,nn))-mat%data_c(nn,n))>tol) call judft_warn("Large unsymmetry in matrix packing")
          endif
          i=i+1
       end DO
    end DO
  end function t_mat_to_packed

  subroutine t_mat_inverse(mat)
    implicit none
    CLASS(t_mat),INTENT(INOUT)       :: mat
    integer                :: info
    real, allocatable      :: work_r(:)
    integer, allocatable   :: ipiv(:)
    complex,allocatable    :: work_c(:)
    
    
    if (mat%matsize1.ne.mat%matsize2) call judft_error("Can only invert square matrices",hint="This is a BUG in FLEUR, please report")
    ALLOCATE(ipiv(mat%matsize1))

    if (mat%l_real) THEN
       ALLOCATE(work_r(mat%matsize1))
       call dgetrf(mat%matsize1,mat%matsize1,mat%data_r,size(mat%data_r,1),ipiv,info)
       if(info.ne.0) call judft_error("Failed to invert matrix: dpotrf failed.")
       call dgetri(mat%matsize1,mat%data_r,size(mat%data_r,1),ipiv,work_r,size(work_r),info)
       if(info.ne.0) call judft_error("Failed to invert matrix: dpotrf failed.")
    else
       ALLOCATE(work_c(mat%matsize1))
       call zgetrf(mat%matsize1,mat%matsize1,mat%data_c,size(mat%data_c,1),ipiv,info)
       if(info.ne.0) call judft_error("Failed to invert matrix: dpotrf failed.")
       call zgetri(mat%matsize1,mat%data_c,size(mat%data_c,1),ipiv,work_c,size(work_c),info)
       if(info.ne.0) call judft_error("Failed to invert matrix: dpotrf failed.")
    end if
  end subroutine t_mat_inverse
     
Daniel Wortmann's avatar
Daniel Wortmann committed
165
end module m_types_rcmat