LCOV - code coverage report
Current view: top level - types - types_mat.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 49 156 31.4 %
Date: 2019-09-08 04:53:50 Functions: 8 16 50.0 %

          Line data    Source code
       1             : MODULE m_types_mat
       2             :   USE m_judft
       3             :   IMPLICIT NONE
       4             :   PRIVATE
       5             :   !<This is the basic type to store and manipulate real/complex rank-2 matrices
       6             :   !!
       7             :   !! In its simple implementation here, it contains a fields for the matrix-size and
       8             :   !! a real and complex array for the data
       9             :   !! This data-type will be overwritten for distributed matrixes by t_mpimat as defined in types_mpimat.F90
      10             :   
      11             :    TYPE :: t_mat
      12             :      LOGICAL :: l_real                     !>Store either real or complex data
      13             :      INTEGER :: matsize1=-1                !> matsize1=size(data_?,1),i.e. no of rows
      14             :      INTEGER :: matsize2=-1                !> matsize2=size(data_?,2),i.e. no of columns
      15             :      REAL,ALLOCATABLE CPP_MANAGED    :: data_r(:,:)
      16             :      COMPLEX,ALLOCATABLE CPP_MANAGED :: data_c(:,:)
      17             :    CONTAINS
      18             :      PROCEDURE        :: alloc => t_mat_alloc                !> allocate the data-arrays
      19             :      PROCEDURE        :: multiply=>t_mat_multiply            !> do a matrix-matrix multiply
      20             :      PROCEDURE        :: transpose=>t_mat_transpose          !> transpose the matrix
      21             :      PROCEDURE        :: from_packed=>t_mat_from_packed      !> initialized from a packed-storage matrix
      22             :      PROCEDURE        :: inverse =>t_mat_inverse             !> invert the matrix
      23             :      PROCEDURE        :: linear_problem => t_mat_lproblem    !> Solve linear equation
      24             :      PROCEDURE        :: to_packed=>t_mat_to_packed          !> convert to packed-storage matrix
      25             :      PROCEDURE        :: clear => t_mat_clear                !> set data arrays to zero
      26             :      PROCEDURE        :: copy => t_mat_copy                  !> copy into another t_mat (overloaded for t_mpimat)
      27             :      PROCEDURE        :: move => t_mat_move                  !> move data into another t_mat (overloaded for t_mpimat)
      28             :      PROCEDURE        :: init_details => t_mat_init
      29             :      PROCEDURE        :: init_template => t_mat_init_template              !> initalize the matrix(overloaded for t_mpimat)
      30             :      GENERIC          :: init => init_details,init_template
      31             :      PROCEDURE        :: free => t_mat_free                  !> dealloc the data (overloaded for t_mpimat)
      32             :      PROCEDURE        :: add_transpose => t_mat_add_transpose!> add the tranpose/Hermitian conjg. without the diagonal (overloaded for t_mpimat)
      33             :    END type t_mat
      34             :    PUBLIC t_mat
      35             :  CONTAINS
      36             : 
      37           0 :    SUBROUTINE t_mat_lproblem(mat,vec)
      38             :      IMPLICIT NONE
      39             :      CLASS(t_mat),INTENT(IN)     :: mat
      40             :      TYPE(t_mat),INTENT(INOUT)   :: vec
      41             : 
      42             :      INTEGER:: lwork,info
      43           0 :      REAL,ALLOCATABLE:: work(:)
      44           0 :      INTEGER,allocatable::ipiv(:)
      45             :     
      46           0 :      IF ((mat%l_real.NEQV.vec%l_real).OR.(mat%matsize1.NE.mat%matsize2).OR.(mat%matsize1.NE.vec%matsize1)) &
      47           0 :           CALL judft_error("Invalid matices in t_mat_lproblem")
      48           0 :      IF (mat%l_real) THEN
      49           0 :         IF (ALL(ABS(mat%data_r-TRANSPOSE(mat%data_r))<1E-8)) THEN
      50             :            !Matrix is symmetric
      51           0 :            CALL DPOSV( 'Upper', mat%matsize1, vec%matsize2, mat%data_r, mat%matsize1, vec%data_r, vec%matsize1, INFO )
      52           0 :            IF (INFO>0) THEN
      53             :               !Matrix was not positive definite
      54           0 :               lwork=-1;ALLOCATE(work(1))
      55           0 :               CALL DSYSV( 'Upper', mat%matsize1, vec%matsize2, mat%data_r, mat%matsize1, IPIV, vec%data_r, vec%matsize1, WORK, LWORK,INFO )
      56           0 :               lwork=INT(work(1))
      57           0 :               DEALLOCATE(work);ALLOCATE(ipiv(mat%matsize1),work(lwork))
      58           0 :               CALL DSYSV( 'Upper', mat%matsize1, vec%matsize2, mat%data_r, mat%matsize1, IPIV, vec%data_r, vec%matsize1, WORK, LWORK,INFO )
      59           0 :               IF (info.NE.0) CALL judft_error("Could not solve linear equation, matrix singular")
      60             :            END IF
      61             :         ELSE
      62           0 :            CALL judft_error("TODO: mode not implemented in t_mat_lproblem")
      63             :         END IF
      64             :      ELSE
      65           0 :         CALL judft_error("TODO: mode not implemented in t_mat_lproblem")
      66             :      ENDIF
      67           0 :    END SUBROUTINE t_mat_lproblem
      68             : 
      69        1520 :    SUBROUTINE t_mat_free(mat)
      70             :      CLASS(t_mat),INTENT(INOUT)::mat
      71        1520 :      IF (ALLOCATED(mat%data_c)) DEALLOCATE(mat%data_c)
      72        1520 :      IF (ALLOCATED(mat%data_r)) DEALLOCATE(mat%data_r)
      73        1520 :    END SUBROUTINE t_mat_free
      74             :    
      75           0 :    SUBROUTINE t_mat_add_transpose(mat,mat1)
      76             :     CLASS(t_mat),INTENT(INOUT)::mat,mat1
      77             :     INTEGER::i,j
      78           0 :     IF ((mat%matsize1.NE.mat1%matsize2).OR. &
      79             :          (mat%matsize2.NE.mat1%matsize1)) &
      80           0 :          CALL judft_error("Matrix sizes missmatch in add_transpose")
      81           0 :     IF (mat%l_real.AND.mat1%l_real) THEN
      82           0 :        DO i=1,mat%matsize2
      83           0 :           DO j=i+1,mat%matsize1
      84           0 :              mat%data_r(j,i)=mat1%data_r(i,j)
      85             :           ENDDO
      86             :        ENDDO
      87           0 :     ELSEIF((.NOT.mat%l_real).AND.(.NOT.mat1%l_real)) THEN
      88           0 :        DO i=1,mat%matsize2
      89           0 :           DO j=i+1,mat%matsize1
      90           0 :              mat%data_c(j,i)=CONJG(mat1%data_c(i,j))
      91             :           ENDDO
      92             :        ENDDO
      93             :     ELSE
      94           0 :        call judft_error("Inconsistency between data types in m_mat")
      95             :     END IF
      96           0 :   END SUBROUTINE t_mat_add_transpose
      97             : 
      98             :   
      99             :    
     100        2460 :    SUBROUTINE t_mat_init(mat,l_real,matsize1,matsize2,mpi_subcom,l_2d,nb_x,nb_y)
     101             :      CLASS(t_mat) :: mat
     102             :      LOGICAL,INTENT(IN),OPTIONAL:: l_real
     103             :      INTEGER,INTENT(IN),OPTIONAL:: matsize1,matsize2
     104             :      INTEGER,INTENT(IN),OPTIONAL:: mpi_subcom,nb_x,nb_y !not needed here, only for allowing overloading this in t_mpimat
     105             :      LOGICAL,INTENT(IN),OPTIONAL:: l_2d                 !not needed here either
     106             : 
     107        2460 :      CALL mat%alloc(l_real,matsize1,matsize2)
     108        2460 :    END SUBROUTINE t_mat_init
     109         608 :    SUBROUTINE t_mat_init_template(mat,templ,global_size1,global_size2)
     110             :      IMPLICIT NONE
     111             :      CLASS(t_mat),INTENT(INOUT) :: mat
     112             :      CLASS(t_mat),INTENT(IN)    :: templ
     113             :      INTEGER,INTENT(IN),OPTIONAL:: global_size1,global_size2
     114             : 
     115         608 :      IF (PRESENT(global_size1).AND.PRESENT(global_size2)) THEN
     116           0 :         IF ((global_size1.NE.templ%matsize1).OR.(global_size2.NE.templ%matsize2)) CALL judft_error("BUG:Invalid change of size in init by template")
     117             :      END IF
     118         608 :      mat%l_real=templ%l_real
     119         608 :      mat%matsize1=templ%matsize1
     120         608 :      mat%matsize2=templ%matsize2
     121         608 :      IF (mat%l_real) THEN
     122         304 :         ALLOCATE(mat%data_r(mat%matsize1,mat%matsize2))
     123         304 :         ALLOCATE(mat%data_c(1,1))
     124       41084 :         mat%data_r=0.0
     125             :      ELSE
     126         304 :         ALLOCATE(mat%data_c(mat%matsize1,mat%matsize2))
     127         304 :         ALLOCATE(mat%data_r(1,1))
     128       58938 :         mat%data_c=0.0
     129             :      END IF
     130         608 :    END SUBROUTINE t_mat_init_template
     131             :      
     132       14586 :   SUBROUTINE t_mat_alloc(mat,l_real,matsize1,matsize2,init)
     133             :     CLASS(t_mat) :: mat
     134             :     LOGICAL,INTENT(IN),OPTIONAL:: l_real
     135             :     INTEGER,INTENT(IN),OPTIONAL:: matsize1,matsize2
     136             :     REAL,INTENT(IN),OPTIONAL   :: init
     137             : 
     138             :     INTEGER:: err
     139       14586 :     IF (present(l_real)) mat%l_real=l_real
     140       14586 :     IF (present(matsize1)) mat%matsize1=matsize1
     141       14586 :     IF (present(matsize2)) mat%matsize2=matsize2
     142             : 
     143       14586 :     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")
     144       14586 :     IF (mat%matsize2<0)  mat%matsize2=mat%matsize1 !Square by default
     145             :     
     146       14586 :     IF (allocated(mat%data_r)) deallocate(mat%data_r)
     147       14586 :     IF (allocated(mat%data_c)) deallocate(mat%data_c)
     148             :     
     149       14586 :     IF (mat%l_real) THEN
     150        2072 :        ALLOCATE(mat%data_r(mat%matsize1,mat%matsize2),STAT=err)
     151        2072 :        ALLOCATE(mat%data_c(0,0))
     152        2072 :        mat%data_r=0.0
     153        2072 :        if (present(init)) mat%data_r=init
     154             :     ELSE
     155       12514 :        ALLOCATE(mat%data_r(0,0))
     156       12514 :        ALLOCATE(mat%data_c(mat%matsize1,mat%matsize2),STAT=err)
     157       12514 :        mat%data_c=0.0
     158       12514 :        IF (PRESENT(init)) mat%data_c=init
     159             :     ENDIF
     160             : 
     161       14586 :     IF (err>0) CALL judft_error("Allocation of memmory failed for mat datatype",hint="You probably run out of memory")
     162       14586 :   END SUBROUTINE t_mat_alloc
     163             : 
     164           0 :   SUBROUTINE t_mat_multiply(mat1,mat2,res)
     165             :     CLASS(t_mat),INTENT(INOUT)        ::mat1
     166             :     CLASS(t_mat),INTENT(IN)           ::mat2
     167             :     CLASS(t_mat),INTENT(OUT),OPTIONAL ::res
     168             : 
     169           0 :     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")
     170             :     
     171           0 :     IF (present(res)) THEN
     172           0 :        call res%alloc(mat1%l_real,mat1%matsize1,mat2%matsize2)
     173           0 :        IF (mat1%l_real) THEN
     174           0 :           res%data_r=matmul(mat1%data_r(:mat1%matsize1,:mat1%matsize2),mat2%data_r(:mat2%matsize1,:mat2%matsize2))
     175             :        ELSE
     176           0 :           res%data_c=matmul(mat1%data_c(:mat1%matsize1,:mat1%matsize2),mat2%data_c(:mat2%matsize1,:mat2%matsize2))
     177             :        ENDIF
     178             :     else
     179           0 :        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")
     180           0 :        if (mat1%l_real) THEN
     181           0 :           mat1%data_r(:mat1%matsize1,:mat1%matsize2)=matmul(mat1%data_r(:mat1%matsize1,:mat1%matsize2),mat2%data_r(:mat2%matsize1,:mat2%matsize2))
     182             :        ELSE
     183           0 :           mat1%data_c(:mat1%matsize1,:mat1%matsize2)=matmul(mat1%data_c(:mat1%matsize1,:mat1%matsize2),mat2%data_c(:mat2%matsize1,:mat2%matsize2))
     184             :        ENDIF
     185             :     end IF
     186           0 :   end SUBROUTINE t_mat_multiply
     187             : 
     188             : 
     189           0 :   SUBROUTINE t_mat_transpose(mat1,res)
     190             :     CLASS(t_mat),INTENT(INOUT)       ::mat1
     191             :     TYPE(t_mat),INTENT(OUT),OPTIONAL ::res
     192             :     
     193           0 :     IF (present(res)) THEN
     194           0 :        call res%alloc(mat1%l_real,mat1%matsize2,mat1%matsize1)
     195           0 :        IF (mat1%l_real) THEN
     196           0 :           res%data_r=transpose(mat1%data_r(:mat1%matsize1,:mat1%matsize2))
     197             :        ELSE
     198           0 :           res%data_c=transpose(mat1%data_c(:mat1%matsize1,:mat1%matsize2))
     199             :        ENDIF
     200             :     else
     201           0 :        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")
     202           0 :        IF (mat1%l_real) THEN
     203           0 :           mat1%data_r(:mat1%matsize1,:mat1%matsize2)=transpose(mat1%data_r(:mat1%matsize1,:mat1%matsize2))
     204             :        ELSE
     205           0 :           mat1%data_c(:mat1%matsize1,:mat1%matsize2)=transpose(mat1%data_c(:mat1%matsize1,:mat1%matsize2))
     206             :        ENDIF
     207             :     end IF
     208           0 :   end SUBROUTINE t_mat_transpose
     209             : 
     210           0 :   SUBROUTINE t_mat_from_packed(mat1,l_real,matsize,packed_r,packed_c)
     211             :     CLASS(t_mat),INTENT(INOUT)       :: mat1
     212             :     INTEGER,INTENT(IN)               :: matsize
     213             :     LOGICAL,INTENT(IN)               :: l_real
     214             :     REAL,INTENT(IN)                  :: packed_r(:)
     215             :     COMPLEX,INTENT(IN)               :: packed_c(:)
     216             : 
     217             :     INTEGER:: n,nn,i
     218           0 :     call mat1%alloc(l_real,matsize,matsize)
     219           0 :     i=1
     220           0 :     DO n=1,matsize
     221           0 :        DO nn=1,n
     222           0 :           if (l_real) THEN
     223           0 :              mat1%data_r(n,nn)=packed_r(i)
     224           0 :              mat1%data_r(nn,n)=packed_r(i)
     225             :           else
     226           0 :              mat1%data_c(n,nn)=conjg(packed_c(i))
     227           0 :              mat1%data_c(nn,n)=packed_c(i)
     228             :           end if
     229           0 :           i=i+1
     230             :        end DO
     231             :     end DO
     232           0 :   end SUBROUTINE t_mat_from_packed
     233             : 
     234           0 :   function t_mat_to_packed(mat)result(packed)
     235             :     CLASS(t_mat),INTENT(IN)       :: mat
     236             :     COMPLEX                       :: packed(mat%matsize1*(mat%matsize1+1)/2)
     237             :     integer :: n,nn,i
     238             :     real,parameter :: tol=1e-5
     239           0 :     if (mat%matsize1.ne.mat%matsize2) call judft_error("Could not pack no-square matrix",hint='This is a BUG, please report')
     240           0 :     i=1
     241           0 :     DO n=1,mat%matsize1
     242           0 :        DO nn=1,n
     243           0 :           if (mat%l_real) THEN
     244           0 :              packed(i)=(mat%data_r(n,nn)+mat%data_r(nn,n))/2.
     245           0 :              if (abs(mat%data_r(n,nn)-mat%data_r(nn,n))>tol) call judft_warn("Large unsymmetry in matrix packing")
     246             :           else
     247           0 :              packed(i)=(conjg(mat%data_c(n,nn))+mat%data_c(nn,n))/2.
     248           0 :              if (abs(conjg(mat%data_c(n,nn))-mat%data_c(nn,n))>tol) call judft_warn("Large unsymmetry in matrix packing")
     249             :           endif
     250           0 :           i=i+1
     251             :        end DO
     252             :     end DO
     253           0 :   end function t_mat_to_packed
     254             : 
     255           0 :   subroutine t_mat_inverse(mat)
     256             :     implicit none
     257             :     CLASS(t_mat),INTENT(INOUT)       :: mat
     258             :     integer                :: info
     259           0 :     real, allocatable      :: work_r(:)
     260           0 :     integer, allocatable   :: ipiv(:)
     261           0 :     complex,allocatable    :: work_c(:)
     262             :     
     263             :     
     264           0 :     if (mat%matsize1.ne.mat%matsize2) call judft_error("Can only invert square matrices",hint="This is a BUG in FLEUR, please report")
     265           0 :     ALLOCATE(ipiv(mat%matsize1))
     266             : 
     267           0 :     if (mat%l_real) THEN
     268           0 :        ALLOCATE(work_r(mat%matsize1))
     269           0 :        call dgetrf(mat%matsize1,mat%matsize1,mat%data_r,size(mat%data_r,1),ipiv,info)
     270           0 :        if(info.ne.0) call judft_error("Failed to invert matrix: dpotrf failed.")
     271           0 :        call dgetri(mat%matsize1,mat%data_r,size(mat%data_r,1),ipiv,work_r,size(work_r),info)
     272           0 :        if(info.ne.0) call judft_error("Failed to invert matrix: dpotrf failed.")
     273             :     else
     274           0 :        ALLOCATE(work_c(mat%matsize1))
     275           0 :        call zgetrf(mat%matsize1,mat%matsize1,mat%data_c,size(mat%data_c,1),ipiv,info)
     276           0 :        if(info.ne.0) call judft_error("Failed to invert matrix: dpotrf failed.")
     277           0 :        call zgetri(mat%matsize1,mat%data_c,size(mat%data_c,1),ipiv,work_c,size(work_c),info)
     278           0 :        if(info.ne.0) call judft_error("Failed to invert matrix: dpotrf failed.")
     279             :     end if
     280           0 :   end subroutine t_mat_inverse
     281             : 
     282         608 :   SUBROUTINE t_mat_move(mat,mat1)
     283             :     IMPLICIT NONE
     284             :     CLASS(t_mat),INTENT(INOUT):: mat
     285             :     CLASS(t_mat),INTENT(INOUT):: mat1
     286             :     !Special case, the full matrix is copied. Then use move alloc
     287         608 :     IF (mat%l_real) THEN
     288         304 :        CALL move_ALLOC(mat1%data_r,mat%data_r)
     289             :     ELSE
     290         304 :        CALL move_ALLOC(mat1%data_c,mat%data_c)
     291             :     END IF
     292         608 :   END SUBROUTINE t_mat_move
     293             :   
     294           0 :   SUBROUTINE t_mat_copy(mat,mat1,n1,n2)
     295             :     IMPLICIT NONE
     296             :     CLASS(t_mat),INTENT(INOUT):: mat
     297             :     CLASS(t_mat),INTENT(IN)   :: mat1
     298             :     INTEGER,INTENT(IN)        :: n1,n2
     299             : 
     300             :     INTEGER:: i1,i2
     301             : 
     302           0 :     i1=mat%matsize1-n1+1  !space available for first dimension
     303           0 :     i2=mat%matsize2-n2+1
     304           0 :     i1=MIN(i1,mat1%matsize1)
     305           0 :     i2=MIN(i2,mat1%matsize2)
     306           0 :     IF (mat%l_real) THEN
     307           0 :        mat%data_r(n1:n1+i1-1,n2:n2+i2-1)=mat1%data_r(:i1,:i2)
     308             :     ELSE
     309           0 :        mat%data_c(n1:n1+i1-1,n2:n2+i2-1)=mat1%data_c(:i1,:i2)
     310             :     END IF
     311             :        
     312           0 :   END SUBROUTINE t_mat_copy
     313             :  
     314        4992 :   SUBROUTINE t_mat_clear(mat)
     315             :     IMPLICIT NONE
     316             :     CLASS(t_mat),INTENT(INOUT):: mat
     317             : 
     318        4992 :     IF (mat%l_real) THEN
     319           0 :        mat%data_r=0.0
     320             :     ELSE
     321        4992 :        mat%data_c=0.0
     322             :     ENDIF
     323        4992 :   END SUBROUTINE t_mat_clear
     324        7736 : END MODULE m_types_mat
     325             : 
     326             : MODULE m_types_rcmat
     327             :   USE m_types_mat
     328             : END MODULE m_types_rcmat

Generated by: LCOV version 1.13