types_mat.F90 14.8 KB
Newer Older
1
MODULE m_types_mat
2 3 4 5 6 7 8 9 10
   USE m_judft
   IMPLICIT NONE
   PRIVATE
   !<This is the basic type to store and manipulate real/complex rank-2 matrices
   !!
   !! In its simple implementation here, it contains a fields for the matrix-size and
   !! a real and complex array for the data
   !! This data-type will be overwritten for distributed matrixes by t_mpimat as defined in types_mpimat.F90

11
   TYPE :: t_mat
12 13 14 15 16
      LOGICAL :: l_real                     !>Store either real or complex data
      INTEGER :: matsize1 = -1                !> matsize1=size(data_?,1),i.e. no of rows
      INTEGER :: matsize2 = -1                !> matsize2=size(data_?,2),i.e. no of columns
      REAL, ALLOCATABLE CPP_MANAGED    :: data_r(:, :)
      COMPLEX, ALLOCATABLE CPP_MANAGED :: data_c(:, :)
17
   CONTAINS
18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
      PROCEDURE        :: alloc => t_mat_alloc                !> allocate the data-arrays
      PROCEDURE        :: multiply => t_mat_multiply            !> do a matrix-matrix multiply
      PROCEDURE        :: transpose => t_mat_transpose          !> transpose the matrix
      PROCEDURE        :: from_packed => t_mat_from_packed      !> initialized from a packed-storage matrix
      PROCEDURE        :: inverse => t_mat_inverse             !> invert the matrix
      PROCEDURE        :: linear_problem => t_mat_lproblem    !> Solve linear equation
      PROCEDURE        :: to_packed => t_mat_to_packed          !> convert to packed-storage matrix
      PROCEDURE        :: clear => t_mat_clear                !> set data arrays to zero
      PROCEDURE        :: copy => t_mat_copy                  !> copy into another t_mat (overloaded for t_mpimat)
      PROCEDURE        :: move => t_mat_move                  !> move data into another t_mat (overloaded for t_mpimat)
      PROCEDURE        :: init_details => t_mat_init
      PROCEDURE        :: init_template => t_mat_init_template              !> initalize the matrix(overloaded for t_mpimat)
      GENERIC          :: init => init_details, init_template
      PROCEDURE        :: free => t_mat_free                  !> dealloc the data (overloaded for t_mpimat)
      PROCEDURE        :: add_transpose => t_mat_add_transpose!> add the tranpose/Hermitian conjg. without the diagonal (overloaded for t_mpimat)
33
!      PROCEDURE        :: save_npy => t_mat_save_npy
34
   END type t_mat
35
   PUBLIC t_mat
36 37 38 39 40 41 42 43 44 45 46
CONTAINS

   SUBROUTINE t_mat_lproblem(mat, vec)
      IMPLICIT NONE
      CLASS(t_mat), INTENT(IN)     :: mat
      TYPE(t_mat), INTENT(INOUT)   :: vec

      INTEGER:: lwork, info
      REAL, ALLOCATABLE:: work(:)
      INTEGER, allocatable::ipiv(:)

47 48
      IF ((mat%l_real .NEQV. vec%l_real) .OR. (mat%matsize1 .NE. mat%matsize2) &
          .OR. (mat%matsize1 .NE. vec%matsize1)) &
49 50 51 52
         CALL judft_error("Invalid matices in t_mat_lproblem")
      IF (mat%l_real) THEN
         IF (ALL(ABS(mat%data_r - TRANSPOSE(mat%data_r)) < 1E-8)) THEN
            !Matrix is symmetric
53 54
            CALL DPOSV('Upper', mat%matsize1, vec%matsize2, mat%data_r, mat%matsize1,&
                                vec%data_r, vec%matsize1, INFO)
55 56 57
            IF (INFO > 0) THEN
               !Matrix was not positive definite
               lwork = -1; ALLOCATE (work(1))
58 59
               CALL DSYSV('Upper', mat%matsize1, vec%matsize2, mat%data_r, mat%matsize1, IPIV, &
                           vec%data_r, vec%matsize1, WORK, LWORK, INFO)
60 61
               lwork = INT(work(1))
               DEALLOCATE (work); ALLOCATE (ipiv(mat%matsize1), work(lwork))
62 63
               CALL DSYSV('Upper', mat%matsize1, vec%matsize2, mat%data_r, mat%matsize1, IPIV, &
                           vec%data_r, vec%matsize1, WORK, LWORK, INFO)
64 65 66 67 68 69 70 71
               IF (info .NE. 0) CALL judft_error("Could not solve linear equation, matrix singular")
            END IF
         ELSE
            CALL judft_error("TODO: mode not implemented in t_mat_lproblem")
         END IF
      ELSE
         CALL judft_error("TODO: mode not implemented in t_mat_lproblem")
      ENDIF
72 73
   END SUBROUTINE t_mat_lproblem

74
   SUBROUTINE t_mat_free(mat)
75 76 77
      CLASS(t_mat), INTENT(INOUT)::mat
      IF (ALLOCATED(mat%data_c)) DEALLOCATE (mat%data_c)
      IF (ALLOCATED(mat%data_r)) DEALLOCATE (mat%data_r)
78
   END SUBROUTINE t_mat_free
79 80 81 82 83 84

   SUBROUTINE t_mat_add_transpose(mat, mat1)
      CLASS(t_mat), INTENT(INOUT)::mat, mat1
      INTEGER::i, j
      IF ((mat%matsize1 .NE. mat1%matsize2) .OR. &
          (mat%matsize2 .NE. mat1%matsize1)) &
85
         CALL judft_error("Matrix sizes missmatch in add_transpose")
86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
      IF (mat%l_real .AND. mat1%l_real) THEN
         DO i = 1, mat%matsize2
            DO j = i + 1, mat%matsize1
               mat%data_r(j, i) = mat1%data_r(i, j)
            ENDDO
         ENDDO
      ELSEIF ((.NOT. mat%l_real) .AND. (.NOT. mat1%l_real)) THEN
         DO i = 1, mat%matsize2
            DO j = i + 1, mat%matsize1
               mat%data_c(j, i) = CONJG(mat1%data_c(i, j))
            ENDDO
         ENDDO
      ELSE
         call judft_error("Inconsistency between data types in m_mat")
      END IF
   END SUBROUTINE t_mat_add_transpose

   SUBROUTINE t_mat_init(mat, l_real, matsize1, matsize2, mpi_subcom, l_2d, nb_x, nb_y)
      CLASS(t_mat) :: mat
      LOGICAL, INTENT(IN), OPTIONAL:: l_real
      INTEGER, INTENT(IN), OPTIONAL:: matsize1, matsize2
      INTEGER, INTENT(IN), OPTIONAL:: mpi_subcom, nb_x, nb_y !not needed here, only for allowing overloading this in t_mpimat
      LOGICAL, INTENT(IN), OPTIONAL:: l_2d                 !not needed here either

      CALL mat%alloc(l_real, matsize1, matsize2)
111
   END SUBROUTINE t_mat_init
112 113 114 115 116
   SUBROUTINE t_mat_init_template(mat, templ, global_size1, global_size2)
      IMPLICIT NONE
      CLASS(t_mat), INTENT(INOUT) :: mat
      CLASS(t_mat), INTENT(IN)    :: templ
      INTEGER, INTENT(IN), OPTIONAL:: global_size1, global_size2
117

118
      IF (PRESENT(global_size1) .AND. PRESENT(global_size2)) THEN
119
         IF ((global_size1 .NE. templ%matsize1) .OR. (global_size2 .NE. templ%matsize2)) CALL judft_error("BUG:Invalid change of size in init by template")
120 121 122 123 124 125 126 127 128 129 130 131 132
      END IF
      mat%l_real = templ%l_real
      mat%matsize1 = templ%matsize1
      mat%matsize2 = templ%matsize2
      IF (mat%l_real) THEN
         ALLOCATE (mat%data_r(mat%matsize1, mat%matsize2))
         ALLOCATE (mat%data_c(1, 1))
         mat%data_r = 0.0
      ELSE
         ALLOCATE (mat%data_c(mat%matsize1, mat%matsize2))
         ALLOCATE (mat%data_r(1, 1))
         mat%data_c = 0.0
      END IF
133
   END SUBROUTINE t_mat_init_template
134 135 136 137 138 139 140 141 142 143 144

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

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

146
      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")
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
      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)
         ALLOCATE (mat%data_c(0, 0))
         mat%data_r = 0.0
         if (present(init)) mat%data_r = init
      ELSE
         ALLOCATE (mat%data_r(0, 0))
         ALLOCATE (mat%data_c(mat%matsize1, mat%matsize2), STAT=err)
         mat%data_c = 0.0
         IF (PRESENT(init)) mat%data_c = init
      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
      CLASS(t_mat), INTENT(IN)           ::mat2
      CLASS(t_mat), INTENT(OUT), OPTIONAL ::res
171

172
      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")
173 174 175 176 177 178 179 180 181

      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
182
         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")
183
         if (mat1%l_real) THEN
184
            mat1%data_r(:mat1%matsize1, :mat1%matsize2) = matmul(mat1%data_r(:mat1%matsize1, :mat1%matsize2), mat2%data_r(:mat2%matsize1, :mat2%matsize2))
185
         ELSE
186
            mat1%data_c(:mat1%matsize1, :mat1%matsize2) = matmul(mat1%data_c(:mat1%matsize1, :mat1%matsize2), mat2%data_c(:mat2%matsize1, :mat2%matsize2))
187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202
         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
203
         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")
204 205 206 207 208 209 210 211 212
         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

   SUBROUTINE t_mat_from_packed(mat1, l_real, matsize, packed_r, packed_c)
213
      use m_judft
214 215 216 217 218 219 220
      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
221
      call timestart("t_mat_from_packed")
222 223 224 225 226 227 228 229 230 231 232 233 234 235
      call mat1%alloc(l_real, matsize, matsize)
      i = 1
      DO n = 1, matsize
         DO nn = 1, n
            if (l_real) THEN
               mat1%data_r(n, nn) = packed_r(i)
               mat1%data_r(nn, n) = packed_r(i)
            else
               mat1%data_c(n, nn) = conjg(packed_c(i))
               mat1%data_c(nn, n) = packed_c(i)
            end if
            i = i + 1
         end DO
      end DO
236
      call timestop("t_mat_from_packed")
237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267
   end SUBROUTINE t_mat_from_packed

   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(:)

268 269 270 271 272
      call timestart("invert matrix")
      if (mat%matsize1 .ne. mat%matsize2) then
         call judft_error("Can only invert square matrices", &
                          hint="This is a BUG in FLEUR, please report")
      endif
273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
      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
288
      call timestop("invert matrix")
289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332
   end subroutine t_mat_inverse

   SUBROUTINE t_mat_move(mat, mat1)
      IMPLICIT NONE
      CLASS(t_mat), INTENT(INOUT):: mat
      CLASS(t_mat), INTENT(INOUT):: mat1
      !Special case, the full matrix is copied. Then use move alloc
      IF (mat%l_real) THEN
         CALL move_ALLOC(mat1%data_r, mat%data_r)
      ELSE
         CALL move_ALLOC(mat1%data_c, mat%data_c)
      END IF
   END SUBROUTINE t_mat_move

   SUBROUTINE t_mat_copy(mat, mat1, n1, n2)
      IMPLICIT NONE
      CLASS(t_mat), INTENT(INOUT):: mat
      CLASS(t_mat), INTENT(IN)   :: mat1
      INTEGER, INTENT(IN)        :: n1, n2

      INTEGER:: i1, i2

      i1 = mat%matsize1 - n1 + 1  !space available for first dimension
      i2 = mat%matsize2 - n2 + 1
      i1 = MIN(i1, mat1%matsize1)
      i2 = MIN(i2, mat1%matsize2)
      IF (mat%l_real) THEN
         mat%data_r(n1:n1 + i1 - 1, n2:n2 + i2 - 1) = mat1%data_r(:i1, :i2)
      ELSE
         mat%data_c(n1:n1 + i1 - 1, n2:n2 + i2 - 1) = mat1%data_c(:i1, :i2)
      END IF

   END SUBROUTINE t_mat_copy

   SUBROUTINE t_mat_clear(mat)
      IMPLICIT NONE
      CLASS(t_mat), INTENT(INOUT):: mat

      IF (mat%l_real) THEN
         mat%data_r = 0.0
      ELSE
         mat%data_c = 0.0
      ENDIF
   END SUBROUTINE t_mat_clear
Matthias Redies's avatar
Matthias Redies committed
333 334 335 336 337 338 339 340

   subroutine t_mat_save_npy(mat, filename)
      use m_judft
      implicit NONE
      class(t_mat), intent(in) :: mat
      character(len=*)         :: filename

      if (mat%l_real) then
341
!         call save_npy(filename, mat%data_r)
Matthias Redies's avatar
Matthias Redies committed
342
      else
343
!         call save_npy(filename, mat%data_c)
Matthias Redies's avatar
Matthias Redies committed
344 345
      endif
   end subroutine t_mat_save_npy
346
END MODULE m_types_mat