chase_diag.F90 19.3 KB
Newer Older
1 2 3 4 5
! 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
6
! Added MPI implementation, DW 2018
7 8 9
!--------------------------------------------------------------------------------
MODULE m_chase_diag
#ifdef CPP_CHASE
10
  USE m_judft
11 12 13 14 15 16 17 18 19 20 21 22
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

23
  interface
24 25 26 27 28 29 30 31 32
    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

33 34
  !MPI
  INTERFACE
35
     SUBROUTINE mpi_dchase_init( mpi_comm, n, nev, nex, xoff,yoff,xlen,ylen,npr,npc,myrow,mycol) BIND( c, name = 'dchase_init' )
36
       USE, INTRINSIC :: iso_c_binding
37
       INTEGER(c_int) :: mpi_comm, n, nev, nex, xoff,yoff,xlen,ylen,myrow,mycol,npr,npc
38 39 40
     END SUBROUTINE mpi_dchase_init
  END INTERFACE

41
  INTERFACE
42
     SUBROUTINE mpi_zchase_init( mpi_comm, n, nev, nex, xoff,yoff,xlen,ylen,npr,npc,myrow,mycol) BIND( c, name = 'zchase_init' )
43
       USE, INTRINSIC :: iso_c_binding
44
       INTEGER(c_int) :: mpi_comm, n, nev, nex, xoff,yoff,xlen,ylen,myrow,mycol,npr,npc
45 46
     END SUBROUTINE mpi_zchase_init
  END INTERFACE
47 48

  INTERFACE
49
     SUBROUTINE mpi_chase_r(h, v, ritzv, deg, tol, mode, opt ) BIND( c, name = 'dchase_solve' )
50 51 52 53 54 55 56 57
       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

58 59

  INTERFACE
60
     SUBROUTINE mpi_chase_c(h, v, ritzv, deg, tol, mode, opt ) BIND( c, name = 'zchase_solve' )
61 62 63 64 65 66 67 68
       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

69 70

  PRIVATE
71

72
  INTEGER         :: chase_eig_id
Daniel Wortmann's avatar
Daniel Wortmann committed
73
  PUBLIC init_chase
74 75 76
#endif
  REAL            :: scale_distance
  REAL            :: tol
77

Daniel Wortmann's avatar
Daniel Wortmann committed
78
  PUBLIC chase_distance,chase_diag
79 80

CONTAINS
81

82 83 84
  SUBROUTINE chase_distance(dist)
    IMPLICIT NONE
    REAL,INTENT(in)::dist
85

86 87
    tol=MAX(1E-8,dist*scale_distance)
  END SUBROUTINE chase_distance
88

89
#ifdef CPP_CHASE
90
    SUBROUTINE init_chase(mpi,input,atoms,kpts,noco,l_real)
91
    USE m_types_mpimat
92
    USE m_types_setup
93 94 95 96 97 98 99
    USE m_types_mpi
    USE m_judft
    USE m_eig66_io

    IMPLICIT NONE

    TYPE(t_mpi),               INTENT(IN)    :: mpi
100

Daniel Wortmann's avatar
Daniel Wortmann committed
101
    TYPE(t_input),             INTENT(IN)    :: input
102 103 104 105 106 107
    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
108
    CHARACTER(len=1000)::arg
109 110


111
    scale_distance=1E-3
112 113 114 115
    !IF (judft_was_argument("-chase_tol_scale")) THEN
    !   arg=juDFT_string_for_argument("-chase_tol_scale")
    !   READ(arg,*) scale_distance
    !ENDIF
116

117
    IF (TRIM(juDFT_string_for_argument("-diag"))=="chase") THEN
118 119 120
       nevd = min(input%neig,lapw%dim_nvd()+atoms%nlotot)
       nexd = min(max(nevd/4, 45),lapw%dim_nvd()+atoms%nlotot-nevd) !dimensioning for workspace
       chase_eig_id=open_eig(mpi%mpi_comm,lapw%dim_nbasfcn(),nevd+nexd,kpts%nkpt,input%jspins,&
121
                             noco%l_noco,.TRUE.,l_real,noco%l_soc,.FALSE.,mpi%n_size)
122 123
    END IF
  END SUBROUTINE init_chase
124
#endif
125

126
   SUBROUTINE chase_diag(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
127 128
    USE m_types_mpimat
    USE m_types_mat
129 130 131 132 133 134 135 136 137 138 139 140 141 142
    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(:)
143
#ifdef CPP_CHASE
144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
    !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
161
#endif
162
  END SUBROUTINE chase_diag
163
#ifdef CPP_CHASE
164
  SUBROUTINE chase_diag_noMPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
165

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

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

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

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

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

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

191 192 193 194 195 196 197 198 199
    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)

200 201 202 203 204 205 206 207 208 209
    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

210
       ! --> now reduce a * z = eig * b * z to the standard form a' * z' = eig * z'
211 212 213 214 215 216 217 218 219
       ! --> 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

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

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

       ne = nev

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

239
       ! --> recover the generalized eigenvectors z by solving z' = l^t * z
240
       CALL dtrtrs('U','N','N',hmat%matsize1,nev,smat%data_r,smat%matsize1,zMatTemp%data_r,zmat%matsize1,info)
241 242 243 244 245 246 247
       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
248
             zmat%data_r(j,i) = zMatTemp%data_r(j,i)
249 250 251 252 253 254 255 256 257 258 259 260 261 262 263
          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

264
       ! --> now reduce a * z = eig * b * z to the standard form a' * z' = eig * z'
265 266 267 268 269 270 271 272 273
       ! --> 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

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

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

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

       ne = nev

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

294
       ! --> recover the generalized eigenvectors z by solving z' = l^t * z
295
       CALL ztrtrs('U','N','N',hmat%matsize1,nev,smat%data_c,smat%matsize1,zMatTemp%data_c,zmat%matsize1,info)
296 297 298 299 300 301 302
       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
303
             zmat%data_c(j,i) = zMatTemp%data_c(j,i)
304 305 306 307 308 309
          END DO
          eig(i) = eigenvalues(i)
       END DO

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

  SUBROUTINE chase_diag_MPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
    use m_types_mpimat
314
    USE m_types_mat
315 316 317
    USE m_judft
    USE iso_c_binding
    USE m_eig66_io
318

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

322
    TYPE(t_mpimat),            INTENT(INOUT)    :: hmat,smat
323 324 325 326 327 328 329
    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(:)

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

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

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

341
    CALL CPU_TIME(t1)
Uliana Alekseeva's avatar
Uliana Alekseeva committed
342 343 344
    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
345

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

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

378
    CALL hmat%generate_full_matrix()
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
    CALL hmat%from_non_dist(zmattemp)
    call zmatTemp%free()
417

418
    ! --> 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
    !Distribute eigenvalues over PEs
439 440 441 442 443 444
    !ne=0
    !DO i=myid+1,nev,np
    !   ne=ne+1
    !   eig(ne)=eigenvalues(i)
    !ENDDO
    eig(:)=eigenvalues(:size(eig))
445 446 447 448 449 450 451 452 453 454

    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

455 456
  END SUBROUTINE chase_diag_MPI

457 458
  SUBROUTINE priv_init_chasempimat(hmat,mat,nev,nex)
    USE m_types_mpimat
459
    USE m_types_mat
460 461 462 463 464 465 466 467 468
    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
469
    INTEGER     :: k,i,j,l
470 471 472 473 474 475 476 477
    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
478
    ALLOCATE(mat%blacsdata)
Uliana Alekseeva's avatar
Uliana Alekseeva committed
479
    mat%blacsdata%mpi_com=hmat%blacsdata%mpi_com
480 481 482
    mat%global_size1=hmat%global_size1
    mat%global_size2=hmat%global_size1
    mat%l_real=hmat%l_real
483

484
    !Determine rank and no of processors
Uliana Alekseeva's avatar
Uliana Alekseeva committed
485 486
    CALL MPI_COMM_RANK(hmat%blacsdata%mpi_com,myid,ierr)
    CALL MPI_COMM_SIZE(hmat%blacsdata%mpi_com,np,ierr)
487 488 489

    !Init ChASE
    IF (mat%l_real) THEN
Uliana Alekseeva's avatar
Uliana Alekseeva committed
490 491
       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)
492
    ELSE
Uliana Alekseeva's avatar
Uliana Alekseeva committed
493 494
       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)
495 496 497
    ENDIF

    !Determine block-sizes
Uliana Alekseeva's avatar
Uliana Alekseeva committed
498 499
    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)
500

501

Uliana Alekseeva's avatar
Uliana Alekseeva committed
502
    ALLOCATE(iusermap(mat%blacsdata%nprow,mat%blacsdata%npcol))
503 504 505
    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
506
    CALL MPI_ALLREDUCE(MPI_IN_PLACE, iusermap, np,MPI_INTEGER,MPI_MAX,mat%blacsdata%mpi_com,ierr)
507
    !Get the Blacs default context
Uliana Alekseeva's avatar
Uliana Alekseeva committed
508
    CALL BLACS_GET(0,0,mat%blacsdata%blacs_desc(2))
509
    ! Create the Grid
Uliana Alekseeva's avatar
Uliana Alekseeva committed
510
    CALL BLACS_GRIDMAP(mat%blacsdata%blacs_desc(2),iusermap,mat%blacsdata%nprow,mat%blacsdata%nprow,mat%blacsdata%npcol)
511

512

513
    !Now create the matrix
Uliana Alekseeva's avatar
Uliana Alekseeva committed
514 515
    mat%matsize1=numroc(mat%global_size1,nbr,myrow,0,mat%blacsdata%nprow)
    mat%matsize2=numroc(mat%global_size1,nbc,mycol,0,mat%blacsdata%npcol)
516 517 518 519 520 521 522 523 524 525 526 527 528
    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
529
    CALL descinit(mat%blacsdata%blacs_desc,mat%global_size1,mat%global_size2,nbr,nbc,0,0,mat%blacsdata%blacs_desc(2),mat%matsize1,ierr)
530 531 532 533 534
    IF (ierr /=0 ) CALL judft_error('Creation of BLACS descriptor failed')

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

535 536


537 538
  END SUBROUTINE priv_init_chasempimat

539

540
#endif
541
  END MODULE m_chase_diag