LCOV - code coverage report
Current view: top level - diagonalization - chase_diag.F90 (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 0 5 0.0 %
Date: 2024-03-28 04:22:06 Functions: 0 2 0.0 %

          Line data    Source code
       1             : ! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       2             : ! This file is part of FLEUR and available as free software under the conditions
       3             : ! of the MIT license as expressed in the LICENSE file in more detail.
       4             : !
       5             : ! @authors: Miriam Hinzen, Gregor Michalicek
       6             : ! Added MPI implementation, DW 2018
       7             : !--------------------------------------------------------------------------------
       8             : MODULE m_chase_diag
       9             : #ifdef CPP_CHASE
      10             :   USE m_judft
      11             :   USE m_constants
      12             : IMPLICIT NONE
      13             : 
      14             :   interface
      15             :     subroutine chase_c( h, n, v, ritzv, nev, nex, deg, tol, mode, opt ) bind( c, name = 'zchase_' )
      16             :       use, intrinsic :: iso_c_binding
      17             :       complex(c_double_complex)     :: h(n,*), v(n,*)
      18             :       integer(c_int)                :: n, deg, nev, nex
      19             :       real(c_double)                :: ritzv(*), tol
      20             :       character(len=1,kind=c_char)  :: mode, opt
      21             :     end subroutine chase_c
      22             :   end interface
      23             : 
      24             :   interface
      25             :     subroutine chase_r( h, n, v, ritzv, nev, nex, deg, tol, mode, opt ) bind( c, name = 'dchase_' )
      26             :       use, intrinsic :: iso_c_binding
      27             :       real(c_double_complex)        :: h(n,*), v(n,*)
      28             :       integer(c_int)                :: n, deg, nev, nex
      29             :       real(c_double)                :: ritzv(*), tol
      30             :       character(len=1,kind=c_char)  :: mode, opt
      31             :     end subroutine chase_r
      32             :   end interface
      33             : 
      34             :   !MPI
      35             :   INTERFACE
      36             :      SUBROUTINE mpi_dchase_init( mpi_comm, n, nev, nex, xoff,yoff,xlen,ylen,npr,npc,myrow,mycol) BIND( c, name = 'dchase_init' )
      37             :        USE, INTRINSIC :: iso_c_binding
      38             :        INTEGER(c_int) :: mpi_comm, n, nev, nex, xoff,yoff,xlen,ylen,myrow,mycol,npr,npc
      39             :      END SUBROUTINE mpi_dchase_init
      40             :   END INTERFACE
      41             : 
      42             :   INTERFACE
      43             :      SUBROUTINE mpi_zchase_init( mpi_comm, n, nev, nex, xoff,yoff,xlen,ylen,npr,npc,myrow,mycol) BIND( c, name = 'zchase_init' )
      44             :        USE, INTRINSIC :: iso_c_binding
      45             :        INTEGER(c_int) :: mpi_comm, n, nev, nex, xoff,yoff,xlen,ylen,myrow,mycol,npr,npc
      46             :      END SUBROUTINE mpi_zchase_init
      47             :   END INTERFACE
      48             : 
      49             :   INTERFACE
      50             :      SUBROUTINE mpi_chase_r(h, v, ritzv, deg, tol, mode, opt ) BIND( c, name = 'dchase_solve' )
      51             :        USE, INTRINSIC :: iso_c_binding
      52             :        REAL(c_double_complex)        :: h(*), v(*)
      53             :        INTEGER(c_int)                :: deg
      54             :        REAL(c_double)                :: ritzv(*), tol
      55             :        CHARACTER(len=1,kind=c_char)  :: mode, opt
      56             :      END SUBROUTINE mpi_chase_r
      57             :   END INTERFACE
      58             : 
      59             : 
      60             :   INTERFACE
      61             :      SUBROUTINE mpi_chase_c(h, v, ritzv, deg, tol, mode, opt ) BIND( c, name = 'zchase_solve' )
      62             :        USE, INTRINSIC :: iso_c_binding
      63             :        COMPLEX(c_double_complex)     :: h(*), v(*)
      64             :        INTEGER(c_int)                :: deg
      65             :        REAL(c_double)                :: ritzv(*), tol
      66             :        CHARACTER(len=1,kind=c_char)  :: mode, opt
      67             :      END SUBROUTINE mpi_chase_c
      68             :   END INTERFACE
      69             : 
      70             : 
      71             :   PRIVATE
      72             : 
      73             :   INTEGER         :: chase_eig_id
      74             :   PUBLIC init_chase
      75             : #endif
      76             :   REAL            :: scale_distance
      77             :   REAL            :: tol
      78             : 
      79             :   PUBLIC chase_distance,chase_diag
      80             : 
      81             : CONTAINS
      82             : 
      83           0 :   SUBROUTINE chase_distance(dist)
      84             :     IMPLICIT NONE
      85             :     REAL,INTENT(in)::dist
      86             : 
      87           0 :     tol=MAX(1E-8,dist*scale_distance)
      88           0 :   END SUBROUTINE chase_distance
      89             : 
      90             : #ifdef CPP_CHASE
      91             :     SUBROUTINE init_chase(mpi,input,atoms,kpts,noco,l_real)
      92             :     USE m_types_mpimat
      93             :     USE m_types_setup
      94             :     USE m_types_mpi
      95             :     USE m_types_lapw
      96             :     USE m_judft
      97             :     USE m_eig66_io
      98             : 
      99             :     IMPLICIT NONE
     100             : 
     101             :     TYPE(t_mpi),               INTENT(IN)    :: mpi
     102             : 
     103             :     TYPE(t_input),             INTENT(IN)    :: input
     104             :     TYPE(t_atoms),             INTENT(IN)    :: atoms
     105             :     TYPE(t_kpts),              INTENT(IN)    :: kpts
     106             :     TYPE(t_noco),              INTENT(IN)    :: noco
     107             :     LOGICAL,                   INTENT(IN)    :: l_real
     108             : 
     109             :     INTEGER             :: nevd, nexd
     110             :     CHARACTER(len=1000) :: arg
     111             :     TYPE(t_lapw)        :: lapw
     112             : 
     113             :     scale_distance=1E-1
     114             :     !IF (judft_was_argument("-chase_tol_scale")) THEN
     115             :     !   arg=juDFT_string_for_argument("-chase_tol_scale")
     116             :     !   READ(arg,*) scale_distance
     117             :     !ENDIF
     118             : 
     119             :     IF (TRIM(juDFT_string_for_argument("-diag"))=="chase") THEN
     120             :        nevd = min(input%neig,lapw%dim_nvd()+atoms%nlotot)
     121             :        nexd = min(max(nevd/4, 45),lapw%dim_nvd()+atoms%nlotot-nevd) !dimensioning for workspace
     122             :        chase_eig_id=open_eig(mpi%mpi_comm,lapw%dim_nbasfcn(),nevd+nexd,kpts%nkpt,input%jspins,&
     123             :                              noco%l_noco,.TRUE.,l_real,noco%l_soc,.FALSE.,mpi%n_size)
     124             :     END IF
     125             :   END SUBROUTINE init_chase
     126             : #endif
     127             : 
     128           0 :    SUBROUTINE chase_diag(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
     129             :     USE m_types_mpimat
     130             :     USE m_types_mat
     131             :     USE m_judft
     132             :     USE iso_c_binding
     133             :     USE m_eig66_io
     134             : 
     135             :     !Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
     136             :     IMPLICIT NONE
     137             : 
     138             :     CLASS(t_mat),              INTENT(INOUT) :: hmat,smat
     139             :     INTEGER,                   INTENT(IN)    :: ikpt
     140             :     INTEGER,                   INTENT(IN)    :: jsp
     141             :     INTEGER,                   INTENT(IN)    :: iter
     142             :     INTEGER,                   INTENT(INOUT) :: ne
     143             :     CLASS(t_mat), ALLOCATABLE, INTENT(OUT)   :: zmat
     144             :     REAL,                      INTENT(OUT)   :: eig(:)
     145             : #ifdef CPP_CHASE
     146             :     !Choose serial or parallel solver
     147             :     SELECT TYPE(hmat)
     148             :     CLASS is (t_mpimat)
     149             :        SELECT TYPE(smat)
     150             :        CLASS is (t_mpimat)
     151             :           CALL chase_diag_MPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
     152             :        CLASS default
     153             :           CALL judft_error("Inconsistent matrix setup")
     154             :        END SELECT
     155             :     CLASS is (t_mat)
     156             :        SELECT TYPE(smat)
     157             :        CLASS is (t_mat)
     158             :           CALL chase_diag_noMPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
     159             :        CLASS default
     160             :           CALL judft_error("Inconsistent matrix setup")
     161             :        END SELECT
     162             :     END SELECT
     163             : #endif
     164           0 :   END SUBROUTINE chase_diag
     165             : #ifdef CPP_CHASE
     166             :   SUBROUTINE chase_diag_noMPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
     167             : 
     168             :     USE m_types_mat
     169             :     USE m_judft
     170             :     USE iso_c_binding
     171             :     USE m_eig66_io
     172             : 
     173             :     !Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
     174             :     IMPLICIT NONE
     175             : 
     176             :     TYPE(t_mat),               INTENT(INOUT) :: hmat,smat
     177             :     INTEGER,                   INTENT(IN)    :: ikpt
     178             :     INTEGER,                   INTENT(IN)    :: jsp
     179             :     INTEGER,                   INTENT(IN)    :: iter
     180             :     INTEGER,                   INTENT(INOUT) :: ne
     181             :     CLASS(t_mat), ALLOCATABLE, INTENT(OUT)   :: zmat
     182             :     REAL,                      INTENT(OUT)   :: eig(:)
     183             : 
     184             :     INTEGER            :: i, j, nev, nex, nbands
     185             :     INTEGER            :: info
     186             : 
     187             :     CLASS(t_Mat),              ALLOCATABLE  :: zMatTemp
     188             :     REAL(c_double),            ALLOCATABLE  :: eigenvalues(:)
     189             : 
     190             :     ALLOCATE(t_mat::zmat)
     191             :     CALL zmat%alloc(hmat%l_real,hmat%matsize1,ne)
     192             : 
     193             :     nev = min(ne,hmat%matsize1)
     194             :     nex = min(max(nev/4, 45), hmat%matsize1-nev) !dimensioning for workspace
     195             : 
     196             :     ALLOCATE(eigenvalues(nev+nex))
     197             :     eigenvalues = 0.0
     198             : 
     199             :     ALLOCATE(t_mat::zmatTemp)
     200             :     CALL zMatTemp%alloc(hmat%l_real,hmat%matsize1,nev+nex)
     201             : 
     202             :     IF (hmat%l_real) THEN
     203             : 
     204             :        ! --> start with Cholesky factorization of b ( so that b = l * l^t)
     205             :        ! --> b is overwritten by l
     206             :        CALL dpotrf('U',smat%matsize1,smat%data_r,SIZE(smat%data_r,1),info)
     207             :        IF (info.NE.0) THEN
     208             :           WRITE (*,*) 'Error in dpotrf: info =',info
     209             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     210             :        ENDIF
     211             : 
     212             :        ! --> now reduce a * z = eig * b * z to the standard form a' * z' = eig * z'
     213             :        ! --> where a' = (l)^-1 * a * (l^t)^-1 and z' = l^t * z
     214             :        CALL dsygst(1,'U',smat%matsize1,hmat%data_r,SIZE(hmat%data_r,1),smat%data_r,SIZE(smat%data_r,1),info)
     215             :        IF (info.NE.0) THEN
     216             :           WRITE (oUnit,*) 'Error in dsygst: info =',info
     217             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     218             :        ENDIF
     219             : 
     220             :        ! --> solve a' * z' = eig * z' for eigenvalues eig between lb und ub
     221             : 
     222             :        zMatTemp%data_r = 0.0
     223             : 
     224             :        do j = 1, hmat%matsize1
     225             :           do i = 1, j
     226             :              hmat%data_r(j,i) = hmat%data_r(i,j)
     227             :           end do
     228             :        end do
     229             :        if(iter.EQ.1) then
     230             :           CALL chase_r(hmat%data_r, hmat%matsize1, zMatTemp%data_r, eigenvalues, nev, nex, 5, scale_distance, 'R', 'S' )
     231             :        else
     232             :           CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
     233             :           CALL chase_r(hmat%data_r, hmat%matsize1, zMatTemp%data_r, eigenvalues, nev, nex, 5, tol, 'A', 'S' )
     234             :        end if
     235             : 
     236             :        ne = nev
     237             : 
     238             :        CALL write_eig(chase_eig_id,ikpt,jsp,nev+nex,nev+nex,&
     239             :             eigenvalues(:(nev+nex)),zmat=zMatTemp)
     240             : 
     241             :        ! --> recover the generalized eigenvectors z by solving z' = l^t * z
     242             :        CALL dtrtrs('U','N','N',hmat%matsize1,nev,smat%data_r,smat%matsize1,zMatTemp%data_r,zmat%matsize1,info)
     243             :        IF (info.NE.0) THEN
     244             :           WRITE (oUnit,*) 'Error in dtrtrs: info =',info
     245             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     246             :        ENDIF
     247             : 
     248             :        DO i = 1, ne
     249             :           DO j = 1, hmat%matsize1
     250             :              zmat%data_r(j,i) = zMatTemp%data_r(j,i)
     251             :           END DO
     252             :           eig(i) = eigenvalues(i)
     253             :        END DO
     254             : 
     255             : 
     256             :     ELSE
     257             : 
     258             :        ! --> start with Cholesky factorization of b ( so that b = l * l^t)
     259             :        ! --> b is overwritten by l
     260             :        CALL zpotrf('U',smat%matsize1,smat%data_c,SIZE(smat%data_c,1),info)
     261             :        IF (info.NE.0) THEN
     262             :           WRITE (*,*) 'Error in zpotrf: info =',info
     263             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     264             :        ENDIF
     265             : 
     266             :        ! --> now reduce a * z = eig * b * z to the standard form a' * z' = eig * z'
     267             :        ! --> where a' = (l)^-1 * a * (l^t)^-1 and z' = l^t * z
     268             :        CALL zhegst(1,'U',smat%matsize1,hmat%data_c,SIZE(hmat%data_c,1),smat%data_c,SIZE(smat%data_c,1),info)
     269             :        IF (info.NE.0) THEN
     270             :           WRITE (oUnit,*) 'Error in zhegst: info =',info
     271             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     272             :        ENDIF
     273             : 
     274             :        ! --> solve a' * z' = eig * z' for eigenvalues eig between lb und ub
     275             : 
     276             :        zMatTemp%data_c = CMPLX(0.0,0.0)
     277             : 
     278             :        do j = 1, hmat%matsize1
     279             :           do i = 1, j
     280             :              hmat%data_c(j,i) = conjg(hmat%data_c(i,j))
     281             :           end do
     282             :        end do
     283             : 
     284             :        if(iter.EQ.1) then
     285             :           CALL chase_c(hmat%data_c, hmat%matsize1, zMatTemp%data_c, eigenvalues, nev, nex, 5, scale_distance, 'R', 'S' )
     286             :        else
     287             :           CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
     288             :           call chase_c(hmat%data_c, hmat%matsize1, zMatTemp%data_c, eigenvalues, nev, nex, 5, tol, 'A', 'S' )
     289             :        end if
     290             : 
     291             :        ne = nev
     292             : 
     293             :        CALL write_eig(chase_eig_id,ikpt,jsp,nev+nex,nev+nex,&
     294             :             eigenvalues(:(nev+nex)),zmat=zMatTemp)
     295             : 
     296             :        ! --> recover the generalized eigenvectors z by solving z' = l^t * z
     297             :        CALL ztrtrs('U','N','N',hmat%matsize1,nev,smat%data_c,smat%matsize1,zMatTemp%data_c,zmat%matsize1,info)
     298             :        IF (info.NE.0) THEN
     299             :           WRITE (oUnit,*) 'Error in ztrtrs: info =',info
     300             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     301             :        ENDIF
     302             : 
     303             :        DO i = 1, ne
     304             :           DO j = 1, hmat%matsize1
     305             :              zmat%data_c(j,i) = zMatTemp%data_c(j,i)
     306             :           END DO
     307             :           eig(i) = eigenvalues(i)
     308             :        END DO
     309             : 
     310             :     ENDIF
     311             :     IF (info.NE.0) CALL judft_error("Diagonalization via ChASE failed", calledby = 'chase_diag')
     312             :   END SUBROUTINE chase_diag_noMPI
     313             : 
     314             :   SUBROUTINE chase_diag_MPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
     315             :     use m_types_mpimat
     316             :     USE m_types_mat
     317             :     USE m_judft
     318             :     USE iso_c_binding
     319             :     USE m_eig66_io
     320             :     USE mpi
     321             : 
     322             :     !Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
     323             :     IMPLICIT NONE
     324             : 
     325             :     TYPE(t_mpimat),            INTENT(INOUT)    :: hmat,smat
     326             :     INTEGER,                   INTENT(IN)       :: ikpt
     327             :     INTEGER,                   INTENT(IN)       :: jsp
     328             :     INTEGER,                   INTENT(IN)       :: iter
     329             :     INTEGER,                   INTENT(INOUT)    :: ne
     330             :     CLASS(t_mat), ALLOCATABLE, INTENT(OUT)      :: zmat
     331             :     REAL,                      INTENT(OUT)      :: eig(:)
     332             : 
     333             :     INTEGER            :: i, j, nev, nex, nbands,xoff,yoff,xlen,ylen,ierr,nb_x,nb_y
     334             :     INTEGER            :: info,myid,np
     335             :     REAL               :: scale !scaling of eigenvalues from scalapack
     336             : 
     337             :     TYPE(t_mat)  :: zMatTemp
     338             :     TYPE(t_mpimat)                              :: chase_mat
     339             :     REAL,                          ALLOCATABLE  :: eigenvalues(:)
     340             : 
     341             :     REAL :: t1,t2,t3,t4
     342             : 
     343             :     CALL CPU_TIME(t1)
     344             :     CALL MPI_COMM_RANK(hmat%blacsdata%mpi_com,myid,info)
     345             :     CALL MPI_COMM_SIZE(hmat%blacsdata%mpi_com,np,info)
     346             :     smat%blacsdata%blacs_desc=hmat%blacsdata%blacs_desc
     347             : 
     348             :     call smat%u2l()
     349             :     call hmat%u2l()
     350             :     !Transform to standard problem using SCALAPACK
     351             :     IF (hmat%l_real) THEN
     352             :        CALL pdpotrf('U',smat%global_size1,smat%data_r,1,1,smat%blacsdata%blacs_desc,info)
     353             :        IF (info.NE.0) THEN
     354             :           WRITE (*,*) 'Error in pdpotrf: info =',info
     355             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     356             :        ENDIF
     357             :         CALL pdsygst(1,'U',smat%global_size1,hmat%data_r,1,1,hmat%blacsdata%blacs_desc,smat%data_r,1,1,smat%blacsdata%blacs_desc,scale,info)
     358             :         IF (ABS(scale-1)>1E-10) call judft_error("Scale parameter not implemented in chase_diag")
     359             :         IF (info.NE.0) THEN
     360             :           WRITE (oUnit,*) 'Error in pdsygst: info =',info
     361             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     362             :        ENDIF
     363             :     ELSE
     364             :        CALL pzpotrf('U',smat%global_size1,smat%data_c,1,1,smat%blacsdata%blacs_desc,info)
     365             :        IF (info.NE.0) THEN
     366             :           WRITE (*,*) 'Error in pzpotrf: info =',info
     367             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     368             :        ENDIF
     369             :         CALL pzhegst(1,'U',smat%global_size1,hmat%data_c,1,1,smat%blacsdata%blacs_desc,smat%data_c,1,1,smat%blacsdata%blacs_desc,scale,info)
     370             :         IF (ABS(scale-1)>1E-10) call judft_error("Scale parameter not implemented in chase_diag")
     371             :         IF (info.NE.0) THEN
     372             :           WRITE (oUnit,*) 'Error in pzhegst: info =',info
     373             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     374             :        ENDIF
     375             :     END IF
     376             : 
     377             :     nev = MIN(ne,hmat%global_size1)
     378             :     nex = min(max(nev/4, 45), hmat%global_size1-nev) !dimensioning for workspace
     379             : 
     380             :     CALL hmat%u2l()
     381             :     CALL priv_init_chasempimat(hmat,chase_mat,nev,nex)
     382             : 
     383             :     !CALL chase_mat%u2l()
     384             :     ALLOCATE(eigenvalues(nev+nex))
     385             :     eigenvalues = 0.0
     386             :     !ALLOCATE(t_mpimat::zmatTemp)
     387             :     CALL zMatTemp%init(hmat%l_real,hmat%global_size1,nev+nex,MPI_COMM_SELF,.TRUE.) !Generate a pseudo-distributed matrix
     388             : 
     389             :     IF (hmat%l_real) THEN
     390             :        IF(iter.EQ.1) THEN
     391             :           CALL CPU_TIME(t2)
     392             :           CALL mpi_chase_r(chase_mat%data_r, zMatTemp%data_r, eigenvalues,  5, 1E-1, 'R', 'S' )
     393             :           CALL CPU_TIME(t3)
     394             :        ELSE
     395             :           CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
     396             :           CALL CPU_TIME(t2)
     397             :           CALL mpi_chase_r(chase_mat%data_r,  zMatTemp%data_r, eigenvalues, 5, tol, 'A', 'S' )
     398             :           CALL CPU_TIME(t3)
     399             :        END IF
     400             :     ELSE
     401             :        IF(iter.EQ.1) THEN
     402             :           CALL CPU_TIME(t2)
     403             :           CALL mpi_chase_c(chase_mat%data_c,  zMatTemp%data_c, eigenvalues,  5, 1E-1, 'R', 'S' )
     404             :           CALL CPU_TIME(t3)
     405             :        ELSE
     406             :           CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
     407             :           CALL CPU_TIME(t2)
     408             :           CALL mpi_chase_c(chase_mat%data_c,  zMatTemp%data_c, eigenvalues,  5, tol, 'A', 'S' )
     409             :           CALL CPU_TIME(t3)
     410             :        END IF
     411             :     ENDIF
     412             :     call chase_mat%free()
     413             :     ne = nev
     414             :     IF (myid==0) CALL write_eig(chase_eig_id,ikpt,jsp,nev+nex,nev+nex,&
     415             :          eigenvalues(:(nev+nex)),zmat=zMatTemp)
     416             : 
     417             :     CALL hmat%from_non_dist(zmattemp)
     418             :     call zmatTemp%free()
     419             : 
     420             :     ! --> recover the generalized eigenvectors z by solving z' = l^t * z
     421             :     IF (smat%l_real) THEN
     422             :        CALL pdtrtrs('U','N','N',hmat%global_size1,hmat%global_size1,smat%data_r,1,1,smat%blacsdata%blacs_desc,&
     423             :             hmat%data_r,1,1,smat%blacsdata%blacs_desc,info)
     424             :     ELSE
     425             :        CALL pztrtrs('U','N','N',hmat%global_size1,hmat%global_size1,smat%data_c,1,1,smat%blacsdata%blacs_desc,&
     426             :             hmat%data_c,1,1,smat%blacsdata%blacs_desc,info)
     427             :     END IF
     428             :     IF (info.NE.0) THEN
     429             :        WRITE (oUnit,*) 'Error in p?trtrs: info =',info
     430             :        CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     431             :     ENDIF
     432             : 
     433             :     !     Redistribute eigvec from ScaLAPACK distribution to each process
     434             :     !
     435             :     ALLOCATE(t_mpimat::zmat)
     436             :     CALL zmat%init(hmat%l_real,hmat%global_size1,hmat%global_size1,hmat%blacsdata%mpi_com,.FALSE.)
     437             :     CALL zmat%copy(hmat,1,1)
     438             : 
     439             :     !Calculate number of EV found, local for PEs
     440             :     ne=0
     441             :     DO i=myid+1,nev,np
     442             :        ne=ne+1
     443             :     !   eig(ne)=eigenvalues(i)
     444             :     ENDDO
     445             :     !Provide eigenvalues for all PEs
     446             :     eig(:)=eigenvalues(:size(eig))
     447             : 
     448             :     CALL CPU_TIME(t4)
     449             : 
     450             :     IF (myid==0) THEN
     451             :        PRINT *,"Chase Prep:",t2-t1
     452             :        PRINT *,"Chase Call:",t3-t2
     453             :        PRINT *,"Chase Post:",t4-t3
     454             :        PRINT *,"Chase Total:",t4-t1
     455             :     ENDIF
     456             : 
     457             :   END SUBROUTINE chase_diag_MPI
     458             : 
     459             :   SUBROUTINE priv_init_chasempimat(hmat,mat,nev,nex)
     460             :     USE m_types_mpimat
     461             :     USE m_types_mat
     462             :     USE mpi
     463             :     IMPLICIT NONE
     464             :     TYPE(t_mpimat),INTENT(INOUT)::hmat,mat
     465             :     INTEGER,INTENT(IN)          :: nev,nex
     466             :     INTEGER::nbc,nbr
     467             : 
     468             :     INTEGER     :: myrow,mycol
     469             :     INTEGER     :: npblacs,np,myid
     470             :     INTEGER     :: rowlen,collen,rowoff,coloff
     471             :     INTEGER     :: k,i,j,l
     472             :     INTEGER     :: ierr
     473             : 
     474             :     INTEGER,ALLOCATABLE :: iblacsnums(:),ihelp(:),iusermap(:,:)
     475             : 
     476             :     EXTERNAL descinit, blacs_get
     477             :     EXTERNAL blacs_pinfo, blacs_gridinit
     478             :     INTEGER,EXTERNAL::numroc,indxl2g
     479             : 
     480             :     ALLOCATE(mat%blacsdata)
     481             :     mat%blacsdata%mpi_com=hmat%blacsdata%mpi_com
     482             :     mat%global_size1=hmat%global_size1
     483             :     mat%global_size2=hmat%global_size1
     484             :     mat%l_real=hmat%l_real
     485             : 
     486             :     !Determine rank and no of processors
     487             :     CALL MPI_COMM_RANK(hmat%blacsdata%mpi_com,myid,ierr)
     488             :     CALL MPI_COMM_SIZE(hmat%blacsdata%mpi_com,np,ierr)
     489             : 
     490             :     !Init ChASE
     491             :     IF (mat%l_real) THEN
     492             :        CALL mpi_dchase_init(hmat%blacsdata%mpi_com,mat%global_size1, nev, nex, rowoff,coloff,rowlen,collen,&
     493             :             mat%blacsdata%nprow,mat%blacsdata%npcol,myrow,mycol)
     494             :     ELSE
     495             :        CALL mpi_zchase_init(hmat%blacsdata%mpi_com,mat%global_size1, nev, nex, rowoff,coloff,rowlen,collen,&
     496             :             mat%blacsdata%nprow,mat%blacsdata%npcol,myrow,mycol)
     497             :     ENDIF
     498             : 
     499             :     !Determine block-sizes
     500             :     CALL MPI_ALLREDUCE(rowlen,nbr,1,MPI_INTEGER,MPI_MAX,mat%blacsdata%mpi_com,ierr)
     501             :     CALL MPI_ALLREDUCE(collen,nbc,1,MPI_INTEGER,MPI_MAX,mat%blacsdata%mpi_com,ierr)
     502             : 
     503             : 
     504             :     ALLOCATE(iusermap(mat%blacsdata%nprow,mat%blacsdata%npcol))
     505             :     iusermap=-2
     506             :     !Get BLACS ranks for all MPI ranks
     507             :     CALL BLACS_PINFO(iusermap(myrow+1,mycol+1),npblacs)  ! iamblacs = local process rank (e.g. myid)
     508             :     CALL MPI_ALLREDUCE(MPI_IN_PLACE, iusermap, np,MPI_INTEGER,MPI_MAX,mat%blacsdata%mpi_com,ierr)
     509             :     !Get the Blacs default context
     510             :     CALL BLACS_GET(0,0,mat%blacsdata%blacs_desc(2))
     511             :     ! Create the Grid
     512             :     CALL BLACS_GRIDMAP(mat%blacsdata%blacs_desc(2),iusermap,mat%blacsdata%nprow,mat%blacsdata%nprow,mat%blacsdata%npcol)
     513             : 
     514             : 
     515             :     !Now create the matrix
     516             :     mat%matsize1=numroc(mat%global_size1,nbr,myrow,0,mat%blacsdata%nprow)
     517             :     mat%matsize2=numroc(mat%global_size1,nbc,mycol,0,mat%blacsdata%npcol)
     518             :     IF (mat%l_real) THEN
     519             :        ALLOCATE(mat%data_r(mat%matsize1,mat%matsize2))
     520             :     ELSE
     521             :        ALLOCATE(mat%data_c(mat%matsize1,mat%matsize2))
     522             :     END IF
     523             :     !Check for consistency
     524             :     IF (mat%matsize1.NE.rowlen.OR.mat%matsize2.NE.collen) THEN
     525             :        PRINT *,myid,"R:",mat%matsize1,rowlen,nbr
     526             :        PRINT *,myid,"C:",mat%matsize2,collen,nbc
     527             :        CALL judft_error("Distribution failed for chase")
     528             :     ENDIF
     529             : 
     530             :     !Create blacs descriptor for chase matrix
     531             :     CALL descinit(mat%blacsdata%blacs_desc,mat%global_size1,mat%global_size2,nbr,nbc,0,0,mat%blacsdata%blacs_desc(2),mat%matsize1,ierr)
     532             :     IF (ierr /=0 ) CALL judft_error('Creation of BLACS descriptor failed')
     533             : 
     534             :     !Copy data from hmat
     535             :     CALL mat%copy(hmat,1,1)
     536             : 
     537             : 
     538             : 
     539             :   END SUBROUTINE priv_init_chasempimat
     540             : 
     541             : 
     542             : #endif
     543             :   END MODULE m_chase_diag

Generated by: LCOV version 1.14