chase_diag.F90 19.4 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!
! @authors: Miriam Hinzen, Gregor Michalicek
7
! Added MPI implementation, DW 2018
8 9 10
!--------------------------------------------------------------------------------
MODULE m_chase_diag
#ifdef CPP_CHASE
11
  USE m_judft
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
IMPLICIT NONE

  interface
    subroutine chase_c( h, n, v, ritzv, nev, nex, deg, tol, mode, opt ) bind( c, name = 'zchase_' )
      use, intrinsic :: iso_c_binding
      complex(c_double_complex)     :: h(n,*), v(n,*)
      integer(c_int)                :: n, deg, nev, nex
      real(c_double)                :: ritzv(*), tol
      character(len=1,kind=c_char)  :: mode, opt
    end subroutine chase_c
  end interface

  interface 
    subroutine chase_r( h, n, v, ritzv, nev, nex, deg, tol, mode, opt ) bind( c, name = 'dchase_' )
      use, intrinsic :: iso_c_binding
      real(c_double_complex)        :: h(n,*), v(n,*)
      integer(c_int)                :: n, deg, nev, nex
      real(c_double)                :: ritzv(*), tol
      character(len=1,kind=c_char)  :: mode, opt
    end subroutine chase_r
  end interface

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 40 41 42
     END SUBROUTINE mpi_dchase_init
  END INTERFACE

  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 47 48 49
     END SUBROUTINE mpi_zchase_init
  END INTERFACE
 
  INTERFACE 
50
     SUBROUTINE mpi_chase_r(h, v, ritzv, deg, tol, mode, opt ) BIND( c, name = 'dchase_solve' )
51 52 53 54 55 56 57 58 59 60
       USE, INTRINSIC :: iso_c_binding
       REAL(c_double_complex)        :: h(*), v(*)
       INTEGER(c_int)                :: deg
       REAL(c_double)                :: ritzv(*), tol
       CHARACTER(len=1,kind=c_char)  :: mode, opt
     END SUBROUTINE mpi_chase_r
  END INTERFACE
 

  INTERFACE 
61
     SUBROUTINE mpi_chase_c(h, v, ritzv, deg, tol, mode, opt ) BIND( c, name = 'zchase_solve' )
62 63 64 65 66 67 68 69 70
       USE, INTRINSIC :: iso_c_binding
       COMPLEX(c_double_complex)     :: h(*), v(*)
       INTEGER(c_int)                :: deg
       REAL(c_double)                :: ritzv(*), tol
       CHARACTER(len=1,kind=c_char)  :: mode, opt
     END SUBROUTINE mpi_chase_c
  END INTERFACE
  

71 72
  PRIVATE 

73
  INTEGER         :: chase_eig_id
Daniel Wortmann's avatar
Daniel Wortmann committed
74
  PUBLIC init_chase
75 76 77 78
#endif
  REAL            :: scale_distance
  REAL            :: tol
  
Daniel Wortmann's avatar
Daniel Wortmann committed
79
  PUBLIC chase_distance,chase_diag
80 81

CONTAINS
82

83 84 85 86 87 88
  SUBROUTINE chase_distance(dist)
    IMPLICIT NONE
    REAL,INTENT(in)::dist
    
    tol=MAX(1E-8,dist*scale_distance)
  END SUBROUTINE chase_distance
89

90
#ifdef CPP_CHASE
Daniel Wortmann's avatar
Daniel Wortmann committed
91
    SUBROUTINE init_chase(mpi,DIMENSION,input,atoms,kpts,noco,l_real)
92
    USE m_types_mpimat
93
    USE m_types_setup
94 95 96 97 98 99 100 101
    USE m_types_mpi
    USE m_judft
    USE m_eig66_io

    IMPLICIT NONE

    TYPE(t_mpi),               INTENT(IN)    :: mpi
    TYPE(t_dimension),         INTENT(IN)    :: dimension
Daniel Wortmann's avatar
Daniel Wortmann committed
102
    TYPE(t_input),             INTENT(IN)    :: input
103 104 105 106 107 108
    TYPE(t_atoms),             INTENT(IN)    :: atoms
    TYPE(t_kpts),              INTENT(IN)    :: kpts
    TYPE(t_noco),              INTENT(IN)    :: noco
    LOGICAL,                   INTENT(IN)    :: l_real

    INTEGER                                  :: nevd, nexd
109 110 111 112
    CHARACTER(len=1000)::arg
    
   
    scale_distance=1E-3
113 114 115 116
    !IF (judft_was_argument("-chase_tol_scale")) THEN
    !   arg=juDFT_string_for_argument("-chase_tol_scale")
    !   READ(arg,*) scale_distance
    !ENDIF
117
    
118
    IF (TRIM(juDFT_string_for_argument("-diag"))=="chase") THEN
119 120
       nevd = min(dimension%neigd,dimension%nvd+atoms%nlotot)
       nexd = min(max(nevd/4, 45),dimension%nvd+atoms%nlotot-nevd) !dimensioning for workspace
121
       chase_eig_id=open_eig(mpi%mpi_comm,DIMENSION%nbasfcn,nevd+nexd,kpts%nkpt,input%jspins,&
122
                             noco%l_noco,.TRUE.,l_real,noco%l_soc,.FALSE.,mpi%n_size)
123 124
    END IF
  END SUBROUTINE init_chase
125 126
#endif
  
127
   SUBROUTINE chase_diag(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
128 129
    USE m_types_mpimat
    USE m_types_mat
130 131 132 133 134 135 136 137 138 139 140 141 142 143
    USE m_judft
    USE iso_c_binding
    USE m_eig66_io

    !Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
    IMPLICIT NONE

    CLASS(t_mat),              INTENT(INOUT) :: hmat,smat
    INTEGER,                   INTENT(IN)    :: ikpt
    INTEGER,                   INTENT(IN)    :: jsp
    INTEGER,                   INTENT(IN)    :: iter
    INTEGER,                   INTENT(INOUT) :: ne
    CLASS(t_mat), ALLOCATABLE, INTENT(OUT)   :: zmat
    REAL,                      INTENT(OUT)   :: eig(:)
144
#ifdef CPP_CHASE
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
    !Choose serial or parallel solver
    SELECT TYPE(hmat)
    CLASS is (t_mpimat)
       SELECT TYPE(smat)
       CLASS is (t_mpimat)
          CALL chase_diag_MPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
       CLASS default
          CALL judft_error("Inconsistent matrix setup")
       END SELECT
    CLASS is (t_mat)
       SELECT TYPE(smat)
       CLASS is (t_mat)
          CALL chase_diag_noMPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
       CLASS default
          CALL judft_error("Inconsistent matrix setup")
       END SELECT
    END SELECT
162
#endif    
163
  END SUBROUTINE chase_diag
164
#ifdef CPP_CHASE  
165
  SUBROUTINE chase_diag_noMPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
166

167
    USE m_types_mat
168 169
    USE m_judft
    USE iso_c_binding
170
    USE m_eig66_io
171 172 173

    !Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
    IMPLICIT NONE
174

175
    TYPE(t_mat),               INTENT(INOUT) :: hmat,smat
176 177
    INTEGER,                   INTENT(IN)    :: ikpt
    INTEGER,                   INTENT(IN)    :: jsp
178
    INTEGER,                   INTENT(IN)    :: iter
179 180 181 182
    INTEGER,                   INTENT(INOUT) :: ne
    CLASS(t_mat), ALLOCATABLE, INTENT(OUT)   :: zmat
    REAL,                      INTENT(OUT)   :: eig(:)

183
    INTEGER            :: i, j, nev, nex, nbands
184 185
    INTEGER            :: info

186
    CLASS(t_Mat),              ALLOCATABLE  :: zMatTemp
187 188 189 190 191
    REAL(c_double),            ALLOCATABLE  :: eigenvalues(:)

    ALLOCATE(t_mat::zmat)
    CALL zmat%alloc(hmat%l_real,hmat%matsize1,ne)

192 193 194 195 196 197 198 199 200
    nev = min(ne,hmat%matsize1)
    nex = min(max(nev/4, 45), hmat%matsize1-nev) !dimensioning for workspace

    ALLOCATE(eigenvalues(nev+nex))
    eigenvalues = 0.0

    ALLOCATE(t_mat::zmatTemp)
    CALL zMatTemp%alloc(hmat%l_real,hmat%matsize1,nev+nex)

201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
    IF (hmat%l_real) THEN

       ! --> start with Cholesky factorization of b ( so that b = l * l^t)
       ! --> b is overwritten by l
       CALL dpotrf('U',smat%matsize1,smat%data_r,SIZE(smat%data_r,1),info)
       IF (info.NE.0) THEN
          WRITE (*,*) 'Error in dpotrf: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF

       ! --> now reduce a * z = eig * b * z to the standard form a' * z' = eig * z' 
       ! --> where a' = (l)^-1 * a * (l^t)^-1 and z' = l^t * z
       CALL dsygst(1,'U',smat%matsize1,hmat%data_r,SIZE(hmat%data_r,1),smat%data_r,SIZE(smat%data_r,1),info)
       IF (info.NE.0) THEN
          WRITE (6,*) 'Error in dsygst: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF

       ! --> solve a' * z' = eig * z' for eigenvalues eig between lb und ub

221
       zMatTemp%data_r = 0.0
222 223 224 225 226 227

       do j = 1, hmat%matsize1
          do i = 1, j
             hmat%data_r(j,i) = hmat%data_r(i,j)
          end do
       end do
228
       if(iter.EQ.1) then
229
          CALL chase_r(hmat%data_r, hmat%matsize1, zMatTemp%data_r, eigenvalues, nev, nex, 25, scale_distance, 'R', 'S' )
230
       else
231
          CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
232
          CALL chase_r(hmat%data_r, hmat%matsize1, zMatTemp%data_r, eigenvalues, nev, nex, 25, tol, 'A', 'S' )
233
       end if
234 235 236

       ne = nev

237
       CALL write_eig(chase_eig_id,ikpt,jsp,nev+nex,nev+nex,&
238
            eigenvalues(:(nev+nex)),zmat=zMatTemp)
239

240
       ! --> recover the generalized eigenvectors z by solving z' = l^t * z
241
       CALL dtrtrs('U','N','N',hmat%matsize1,nev,smat%data_r,smat%matsize1,zMatTemp%data_r,zmat%matsize1,info)
242 243 244 245 246 247 248
       IF (info.NE.0) THEN
          WRITE (6,*) 'Error in dtrtrs: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF

       DO i = 1, ne
          DO j = 1, hmat%matsize1
249
             zmat%data_r(j,i) = zMatTemp%data_r(j,i)
250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274
          END DO
          eig(i) = eigenvalues(i)
       END DO


    ELSE

       ! --> start with Cholesky factorization of b ( so that b = l * l^t)
       ! --> b is overwritten by l
       CALL zpotrf('U',smat%matsize1,smat%data_c,SIZE(smat%data_c,1),info)
       IF (info.NE.0) THEN
          WRITE (*,*) 'Error in zpotrf: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF

       ! --> now reduce a * z = eig * b * z to the standard form a' * z' = eig * z' 
       ! --> where a' = (l)^-1 * a * (l^t)^-1 and z' = l^t * z
       CALL zhegst(1,'U',smat%matsize1,hmat%data_c,SIZE(hmat%data_c,1),smat%data_c,SIZE(smat%data_c,1),info)
       IF (info.NE.0) THEN
          WRITE (6,*) 'Error in zhegst: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF

       ! --> solve a' * z' = eig * z' for eigenvalues eig between lb und ub

275
       zMatTemp%data_c = CMPLX(0.0,0.0)
276 277 278 279 280 281 282

       do j = 1, hmat%matsize1
          do i = 1, j
             hmat%data_c(j,i) = conjg(hmat%data_c(i,j))
          end do
       end do

283
       if(iter.EQ.1) then
284
          CALL chase_c(hmat%data_c, hmat%matsize1, zMatTemp%data_c, eigenvalues, nev, nex, 25, scale_distance, 'R', 'S' )
285
       else
286
          CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
287
          call chase_c(hmat%data_c, hmat%matsize1, zMatTemp%data_c, eigenvalues, nev, nex, 25, tol, 'A', 'S' )
288
       end if
289 290 291

       ne = nev

292
       CALL write_eig(chase_eig_id,ikpt,jsp,nev+nex,nev+nex,&
293
            eigenvalues(:(nev+nex)),zmat=zMatTemp)
294

295
       ! --> recover the generalized eigenvectors z by solving z' = l^t * z
296
       CALL ztrtrs('U','N','N',hmat%matsize1,nev,smat%data_c,smat%matsize1,zMatTemp%data_c,zmat%matsize1,info)
297 298 299 300 301 302 303
       IF (info.NE.0) THEN
          WRITE (6,*) 'Error in ztrtrs: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF

       DO i = 1, ne
          DO j = 1, hmat%matsize1
304
             zmat%data_c(j,i) = zMatTemp%data_c(j,i)
305 306 307 308 309 310
          END DO
          eig(i) = eigenvalues(i)
       END DO

    ENDIF
    IF (info.NE.0) CALL judft_error("Diagonalization via ChASE failed", calledby = 'chase_diag')
311 312 313 314
  END SUBROUTINE chase_diag_noMPI

  SUBROUTINE chase_diag_MPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
    use m_types_mpimat
315
    USE m_types_mat
316 317 318 319 320 321 322
    USE m_judft
    USE iso_c_binding
    USE m_eig66_io
    
    !Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
    IMPLICIT NONE

323
    TYPE(t_mpimat),            INTENT(INOUT)    :: hmat,smat
324 325 326 327 328 329 330
    INTEGER,                   INTENT(IN)       :: ikpt
    INTEGER,                   INTENT(IN)       :: jsp
    INTEGER,                   INTENT(IN)       :: iter
    INTEGER,                   INTENT(INOUT)    :: ne
    CLASS(t_mat), ALLOCATABLE, INTENT(OUT)      :: zmat
    REAL,                      INTENT(OUT)      :: eig(:)

331
    INTEGER            :: i, j, nev, nex, nbands,xoff,yoff,xlen,ylen,ierr,nb_x,nb_y
332
    INTEGER            :: info,myid,np
333
    REAL               :: scale !scaling of eigenvalues from scalapack
334

335
    TYPE(t_mat)  :: zMatTemp
336 337
    TYPE(t_mpimat)                              :: chase_mat
    REAL,                          ALLOCATABLE  :: eigenvalues(:)
338 339

    REAL :: t1,t2,t3,t4
340 341
    include 'mpif.h'

342
    CALL CPU_TIME(t1)
Uliana Alekseeva's avatar
Uliana Alekseeva committed
343 344 345
    CALL MPI_COMM_RANK(hmat%blacsdata%mpi_com,myid,info)
    CALL MPI_COMM_SIZE(hmat%blacsdata%mpi_com,np,info)
    smat%blacsdata%blacs_desc=hmat%blacsdata%blacs_desc
346

Uliana Alekseeva's avatar
Uliana Alekseeva committed
347 348
    call smat%generate_full_matrix()
    call hmat%generate_full_matrix()
349 350
    !Transform to standard problem using SCALAPACK
    IF (hmat%l_real) THEN
Uliana Alekseeva's avatar
Uliana Alekseeva committed
351
       CALL pdpotrf('U',smat%global_size1,smat%data_r,1,1,smat%blacsdata%blacs_desc,info)
352 353 354 355
       IF (info.NE.0) THEN
          WRITE (*,*) 'Error in pdpotrf: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF
Uliana Alekseeva's avatar
Uliana Alekseeva committed
356
        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)
357
        IF (ABS(scale-1)>1E-10) call judft_error("Scale parameter not implemented in chase_diag")
358 359 360 361 362
        IF (info.NE.0) THEN
          WRITE (6,*) 'Error in pdsygst: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF
    ELSE
Uliana Alekseeva's avatar
Uliana Alekseeva committed
363
       CALL pzpotrf('U',smat%global_size1,smat%data_c,1,1,smat%blacsdata%blacs_desc,info)
364 365 366 367
       IF (info.NE.0) THEN
          WRITE (*,*) 'Error in pzpotrf: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF
Uliana Alekseeva's avatar
Uliana Alekseeva committed
368
        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)
369
        IF (ABS(scale-1)>1E-10) call judft_error("Scale parameter not implemented in chase_diag")
370 371 372 373 374
        IF (info.NE.0) THEN
          WRITE (6,*) 'Error in pzhegst: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF
    END IF
375

376
    nev = MIN(ne,hmat%global_size1)
377
    nex = min(max(nev/4, 45), hmat%global_size1-nev) !dimensioning for workspace
378 379

    CALL priv_init_chasempimat(hmat,chase_mat,nev,nex)
380
    
381
    CALL chase_mat%generate_full_matrix()
382 383
    ALLOCATE(eigenvalues(nev+nex))
    eigenvalues = 0.0
384
    !ALLOCATE(t_mpimat::zmatTemp)
385
    CALL zMatTemp%init(hmat%l_real,hmat%global_size1,nev+nex,MPI_COMM_SELF,.TRUE.) !Generate a pseudo-distributed matrix
386
    
387 388
    IF (hmat%l_real) THEN
       IF(iter.EQ.1) THEN
389 390 391
          CALL CPU_TIME(t2)
          CALL mpi_chase_r(chase_mat%data_r, zMatTemp%data_r, eigenvalues,  25, 1E-4, 'R', 'S' )
          CALL CPU_TIME(t3)
392 393
       ELSE
          CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
394 395 396
          CALL CPU_TIME(t2)
          CALL mpi_chase_r(chase_mat%data_r,  zMatTemp%data_r, eigenvalues, 25, tol, 'A', 'S' )
          CALL CPU_TIME(t3)
397 398 399
       END IF
    ELSE
       IF(iter.EQ.1) THEN
400 401 402
          CALL CPU_TIME(t2)
          CALL mpi_chase_c(chase_mat%data_c,  zMatTemp%data_c, eigenvalues,  25, 1E-4, 'R', 'S' )
          CALL CPU_TIME(t3)
403 404
       ELSE
          CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
405 406 407
          CALL CPU_TIME(t2)
          CALL mpi_chase_c(chase_mat%data_c,  zMatTemp%data_c, eigenvalues,  25, tol, 'A', 'S' )
          CALL CPU_TIME(t3)
408 409
       END IF
    ENDIF
410
    
411
    ne = nev
412 413
    IF (myid==0) CALL write_eig(chase_eig_id,ikpt,jsp,nev+nex,nev+nex,&
         eigenvalues(:(nev+nex)),zmat=zMatTemp)
414

415 416 417 418
    CALL hmat%from_non_dist(zmattemp)
    call zmatTemp%free()
    
    ! --> recover the generalized eigenvectors z by solving z' = l^t * z
419
    IF (smat%l_real) THEN
Uliana Alekseeva's avatar
Uliana Alekseeva committed
420 421
       CALL pdtrtrs('U','N','N',hmat%global_size1,hmat%global_size1,smat%data_r,1,1,smat%blacsdata%blacs_desc,&
            hmat%data_r,1,1,smat%blacsdata%blacs_desc,info)
422
    ELSE
Uliana Alekseeva's avatar
Uliana Alekseeva committed
423 424
       CALL pztrtrs('U','N','N',hmat%global_size1,hmat%global_size1,smat%data_c,1,1,smat%blacsdata%blacs_desc,&
            hmat%data_c,1,1,smat%blacsdata%blacs_desc,info)
425 426 427 428
    END IF
    IF (info.NE.0) THEN
       WRITE (6,*) 'Error in p?trtrs: info =',info
       CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
429 430 431 432 433 434
    ENDIF

    !     Redistribute eigvec from ScaLAPACK distribution to each process
    !     having all eigenvectors corresponding to his eigenvalues as above
    !
    ALLOCATE(t_mpimat::zmat)
Uliana Alekseeva's avatar
Uliana Alekseeva committed
435
    CALL zmat%init(hmat%l_real,hmat%global_size1,hmat%global_size1,hmat%blacsdata%mpi_com,.FALSE.)
436 437
    CALL zmat%copy(hmat,1,1)

438 439 440 441 442 443
    !Distribute eigenvalues over PEs
    ne=0
    DO i=myid+1,nev,np
       ne=ne+1
       eig(ne)=eigenvalues(i)
    ENDDO
444 445 446 447 448 449 450 451 452 453

    CALL CPU_TIME(t4)

    IF (myid==0) THEN
       PRINT *,"Chase Prep:",t2-t1
       PRINT *,"Chase Call:",t3-t2
       PRINT *,"Chase Post:",t4-t3
       PRINT *,"Chase Total:",t4-t1
    ENDIF

454 455
  END SUBROUTINE chase_diag_MPI

456 457
  SUBROUTINE priv_init_chasempimat(hmat,mat,nev,nex)
    USE m_types_mpimat
458
    USE m_types_mat
459 460 461 462 463 464 465 466 467
    IMPLICIT NONE
    TYPE(t_mpimat),INTENT(INOUT)::hmat,mat
    INTEGER,INTENT(IN)          :: nev,nex
    INTEGER::nbc,nbr

    INCLUDE 'mpif.h'
    INTEGER     :: myrow,mycol
    INTEGER     :: npblacs,np,myid
    INTEGER     :: rowlen,collen,rowoff,coloff
468
    INTEGER     :: k,i,j,l
469 470 471 472 473 474 475 476
    INTEGER     :: ierr

    INTEGER,ALLOCATABLE :: iblacsnums(:),ihelp(:),iusermap(:,:)

    EXTERNAL descinit, blacs_get
    EXTERNAL blacs_pinfo, blacs_gridinit
    INTEGER,EXTERNAL::numroc,indxl2g

Uliana Alekseeva's avatar
Uliana Alekseeva committed
477
    ALLOCATE(mat%blacsdata)
Uliana Alekseeva's avatar
Uliana Alekseeva committed
478
    mat%blacsdata%mpi_com=hmat%blacsdata%mpi_com
479 480 481 482 483
    mat%global_size1=hmat%global_size1
    mat%global_size2=hmat%global_size1
    mat%l_real=hmat%l_real
   
    !Determine rank and no of processors
Uliana Alekseeva's avatar
Uliana Alekseeva committed
484 485
    CALL MPI_COMM_RANK(hmat%blacsdata%mpi_com,myid,ierr)
    CALL MPI_COMM_SIZE(hmat%blacsdata%mpi_com,np,ierr)
486 487 488

    !Init ChASE
    IF (mat%l_real) THEN
Uliana Alekseeva's avatar
Uliana Alekseeva committed
489 490
       CALL mpi_dchase_init(hmat%blacsdata%mpi_com,mat%global_size1, nev, nex, rowoff,coloff,rowlen,collen,&
            mat%blacsdata%nprow,mat%blacsdata%npcol,myrow,mycol)
491
    ELSE
Uliana Alekseeva's avatar
Uliana Alekseeva committed
492 493
       CALL mpi_zchase_init(hmat%blacsdata%mpi_com,mat%global_size1, nev, nex, rowoff,coloff,rowlen,collen,&
            mat%blacsdata%nprow,mat%blacsdata%npcol,myrow,mycol)
494 495 496
    ENDIF

    !Determine block-sizes
Uliana Alekseeva's avatar
Uliana Alekseeva committed
497 498
    CALL MPI_ALLREDUCE(rowlen,nbr,1,MPI_INTEGER,MPI_MAX,mat%blacsdata%mpi_com,ierr)
    CALL MPI_ALLREDUCE(collen,nbc,1,MPI_INTEGER,MPI_MAX,mat%blacsdata%mpi_com,ierr)
499 500

    
Uliana Alekseeva's avatar
Uliana Alekseeva committed
501
    ALLOCATE(iusermap(mat%blacsdata%nprow,mat%blacsdata%npcol))
502 503 504
    iusermap=-2
    !Get BLACS ranks for all MPI ranks
    CALL BLACS_PINFO(iusermap(myrow+1,mycol+1),npblacs)  ! iamblacs = local process rank (e.g. myid)
Uliana Alekseeva's avatar
Uliana Alekseeva committed
505
    CALL MPI_ALLREDUCE(MPI_IN_PLACE, iusermap, np,MPI_INTEGER,MPI_MAX,mat%blacsdata%mpi_com,ierr)
506
    !Get the Blacs default context
Uliana Alekseeva's avatar
Uliana Alekseeva committed
507
    CALL BLACS_GET(0,0,mat%blacsdata%blacs_desc(2))
508
    ! Create the Grid
Uliana Alekseeva's avatar
Uliana Alekseeva committed
509
    CALL BLACS_GRIDMAP(mat%blacsdata%blacs_desc(2),iusermap,mat%blacsdata%nprow,mat%blacsdata%nprow,mat%blacsdata%npcol)
510

511
   
512
    !Now create the matrix
Uliana Alekseeva's avatar
Uliana Alekseeva committed
513 514
    mat%matsize1=numroc(mat%global_size1,nbr,myrow,0,mat%blacsdata%nprow)
    mat%matsize2=numroc(mat%global_size1,nbc,mycol,0,mat%blacsdata%npcol)
515 516 517 518 519 520 521 522 523 524 525 526 527
    IF (mat%l_real) THEN
       ALLOCATE(mat%data_r(mat%matsize1,mat%matsize2))
    ELSE
       ALLOCATE(mat%data_c(mat%matsize1,mat%matsize2))
    END IF
    !Check for consistency
    IF (mat%matsize1.NE.rowlen.OR.mat%matsize2.NE.collen) THEN
       PRINT *,myid,"R:",mat%matsize1,rowlen,nbr
       PRINT *,myid,"C:",mat%matsize2,collen,nbc
       CALL judft_error("Distribution failed for chase")
    ENDIF

    !Create blacs descriptor for chase matrix
Uliana Alekseeva's avatar
Uliana Alekseeva committed
528
    CALL descinit(mat%blacsdata%blacs_desc,mat%global_size1,mat%global_size2,nbr,nbc,0,0,mat%blacsdata%blacs_desc(2),mat%matsize1,ierr)
529 530 531 532 533 534 535 536 537
    IF (ierr /=0 ) CALL judft_error('Creation of BLACS descriptor failed')

    !Copy data from hmat
    CALL mat%copy(hmat,1,1)

    
    
  END SUBROUTINE priv_init_chasempimat

538

539
#endif
540
  END MODULE m_chase_diag