chase_diag.F90 18.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 73 74 75 76
  PRIVATE 

    INTEGER :: chase_eig_id

  PUBLIC init_chase, chase_diag

77 78
  CONTAINS

79
  SUBROUTINE init_chase(mpi,dimension,atoms,kpts,noco,l_real)
80
    USE m_types_mpimat
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
    USE m_types
    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
    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

    IF (juDFT_was_argument("-diag:chase")) THEN
       nevd = min(dimension%neigd,dimension%nvd+atoms%nlotot)
       nexd = min(max(nevd/4, 45),dimension%nvd+atoms%nlotot-nevd) !dimensioning for workspace
101 102
       chase_eig_id=open_eig(mpi%mpi_comm,DIMENSION%nbasfcn,nevd+nexd,kpts%nkpt,DIMENSION%jspd,&
                             noco%l_noco,.TRUE.,l_real,noco%l_soc,.FALSE.,mpi%n_size)
103 104 105
    END IF
  END SUBROUTINE init_chase

106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
   SUBROUTINE chase_diag(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
     USE m_types_mpimat
    USE m_types
    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(:)

    !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
  END SUBROUTINE chase_diag
  
Daniel Wortmann's avatar
Daniel Wortmann committed
143
  SUBROUTINE chase_diag_noMPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
144 145 146 147

    USE m_types
    USE m_judft
    USE iso_c_binding
148
    USE m_eig66_io
149 150 151

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

153
    TYPE(t_mat),               INTENT(INOUT) :: hmat,smat
154 155
    INTEGER,                   INTENT(IN)    :: ikpt
    INTEGER,                   INTENT(IN)    :: jsp
156
    INTEGER,                   INTENT(IN)    :: iter
157 158 159 160
    INTEGER,                   INTENT(INOUT) :: ne
    CLASS(t_mat), ALLOCATABLE, INTENT(OUT)   :: zmat
    REAL,                      INTENT(OUT)   :: eig(:)

161
    INTEGER            :: i, j, nev, nex, nbands
162 163
    INTEGER            :: info

164
    CLASS(t_Mat),              ALLOCATABLE  :: zMatTemp
165 166 167 168 169
    REAL(c_double),            ALLOCATABLE  :: eigenvalues(:)

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

170 171 172 173 174 175 176 177 178
    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)

179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
    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

199
       zMatTemp%data_r = 0.0
200 201 202 203 204 205

       do j = 1, hmat%matsize1
          do i = 1, j
             hmat%data_r(j,i) = hmat%data_r(i,j)
          end do
       end do
206
       if(iter.EQ.1) then
207
          call chase_r(hmat%data_r, hmat%matsize1, zMatTemp%data_r, eigenvalues, nev, nex, 25, 1e-6, 'R', 'S' )
208
       else
209
          CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
210
          call chase_r(hmat%data_r, hmat%matsize1, zMatTemp%data_r, eigenvalues, nev, nex, 25, 1e-6, 'A', 'S' )
211
       end if
212 213 214

       ne = nev

215
       CALL write_eig(chase_eig_id,ikpt,jsp,nev+nex,nev+nex,&
Daniel Wortmann's avatar
Daniel Wortmann committed
216
            eigenvalues(:(nev+nex)),zmat=zMatTemp)
217

218
       ! --> recover the generalized eigenvectors z by solving z' = l^t * z
219
       CALL dtrtrs('U','N','N',hmat%matsize1,nev,smat%data_r,smat%matsize1,zMatTemp%data_r,zmat%matsize1,info)
220 221 222 223 224 225 226
       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
227
             zmat%data_r(j,i) = zMatTemp%data_r(j,i)
228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252
          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

253
       zMatTemp%data_c = CMPLX(0.0,0.0)
254 255 256 257 258 259 260

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

261
       if(iter.EQ.1) then
262
          call chase_c(hmat%data_c, hmat%matsize1, zMatTemp%data_c, eigenvalues, nev, nex, 25, 1e-6, 'R', 'S' )
263
       else
264
          CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
265
          call chase_c(hmat%data_c, hmat%matsize1, zMatTemp%data_c, eigenvalues, nev, nex, 25, 1e-6, 'A', 'S' )
266
       end if
267 268 269

       ne = nev

270
       CALL write_eig(chase_eig_id,ikpt,jsp,nev+nex,nev+nex,&
Daniel Wortmann's avatar
Daniel Wortmann committed
271
            eigenvalues(:(nev+nex)),zmat=zMatTemp)
272

273
       ! --> recover the generalized eigenvectors z by solving z' = l^t * z
274
       CALL ztrtrs('U','N','N',hmat%matsize1,nev,smat%data_c,smat%matsize1,zMatTemp%data_c,zmat%matsize1,info)
275 276 277 278 279 280 281
       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
282
             zmat%data_c(j,i) = zMatTemp%data_c(j,i)
283 284 285 286 287 288
          END DO
          eig(i) = eigenvalues(i)
       END DO

    ENDIF
    IF (info.NE.0) CALL judft_error("Diagonalization via ChASE failed", calledby = 'chase_diag')
289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308
  END SUBROUTINE chase_diag_noMPI

  SUBROUTINE chase_diag_MPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
    use m_types_mpimat
    USE m_types
    USE m_judft
    USE iso_c_binding
    USE m_eig66_io
    
    !Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
    IMPLICIT NONE

    TYPE(t_mpimat),               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(:)

309
    INTEGER            :: i, j, nev, nex, nbands,xoff,yoff,xlen,ylen,ierr,nb_x,nb_y
Daniel Wortmann's avatar
Daniel Wortmann committed
310
    INTEGER            :: info,myid,np
311
    REAL               :: scale !scaling of eigenvalues from scalapack
312

313 314 315
    CLASS(t_mat),                  ALLOCATABLE  :: zMatTemp
    TYPE(t_mpimat)                              :: chase_mat
    REAL,                          ALLOCATABLE  :: eigenvalues(:)
316 317
    include 'mpif.h'

Daniel Wortmann's avatar
Daniel Wortmann committed
318 319 320 321 322 323
    CALL MPI_COMM_RANK(hmat%mpi_com,myid,info)
    CALL MPI_COMM_SIZE(hmat%mpi_com,np,info)
  smat%blacs_desc=hmat%blacs_desc

  call smat%generate_full_matrix()
  call hmat%generate_full_matrix()
324 325 326 327 328 329 330
    !Transform to standard problem using SCALAPACK
    IF (hmat%l_real) THEN
       CALL pdpotrf('U',smat%global_size1,smat%data_r,1,1,smat%blacs_desc,info)
       IF (info.NE.0) THEN
          WRITE (*,*) 'Error in pdpotrf: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
331
        CALL pdsygst(1,'U',smat%global_size1,hmat%data_r,1,1,hmat%blacs_desc,smat%data_r,1,1,smat%blacs_desc,scale,info)
332
        IF (ABS(scale-1)>1E-10) call judft_error("Scale parameter not implemented in chase_diag")
333 334 335 336 337 338 339 340 341 342
        IF (info.NE.0) THEN
          WRITE (6,*) 'Error in pdsygst: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF
    ELSE
       CALL pzpotrf('U',smat%global_size1,smat%data_c,1,1,smat%blacs_desc,info)
       IF (info.NE.0) THEN
          WRITE (*,*) 'Error in pzpotrf: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF
343 344
        CALL pzhegst(1,'U',smat%global_size1,hmat%data_c,1,1,smat%blacs_desc,smat%data_c,1,1,smat%blacs_desc,scale,info)
        IF (ABS(scale-1)>1E-10) call judft_error("Scale parameter not implemented in chase_diag")
345 346 347 348 349
        IF (info.NE.0) THEN
          WRITE (6,*) 'Error in pzhegst: info =',info
          CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
       ENDIF
    END IF
350

Daniel Wortmann's avatar
Daniel Wortmann committed
351
    nev = MIN(ne,hmat%global_size1)
352
    nex = min(max(nev/4, 45), hmat%global_size1-nev) !dimensioning for workspace
353 354

    CALL priv_init_chasempimat(hmat,chase_mat,nev,nex)
355
    
Daniel Wortmann's avatar
Daniel Wortmann committed
356
    CALL chase_mat%generate_full_matrix()
357 358 359
    ALLOCATE(eigenvalues(nev+nex))
    eigenvalues = 0.0
    ALLOCATE(t_mpimat::zmatTemp)
360
    CALL zMatTemp%init(hmat%l_real,hmat%global_size1,nev+nex,MPI_COMM_SELF,.TRUE.) !Generate a pseudo-distributed matrix
361 362 363

    IF (hmat%l_real) THEN
       IF(iter.EQ.1) THEN
364
          CALL mpi_chase_r(chase_mat%data_r, zMatTemp%data_r, eigenvalues,  25, 1e-10, 'R', 'S' )
365 366
       ELSE
          CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
367
          CALL mpi_chase_r(chase_mat%data_r,  zMatTemp%data_r, eigenvalues, 25, 1e-10, 'A', 'S' )
368 369 370
       END IF
    ELSE
       IF(iter.EQ.1) THEN
371
          CALL mpi_chase_c(chase_mat%data_c,  zMatTemp%data_c, eigenvalues,  25, 1e-10, 'R', 'S' )
372 373
       ELSE
          CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
374
          CALL mpi_chase_c(chase_mat%data_c,  zMatTemp%data_c, eigenvalues,  25, 1e-10, 'A', 'S' )
375 376
       END IF
    ENDIF
377

378
    ne = nev
Daniel Wortmann's avatar
Daniel Wortmann committed
379 380
    IF (myid==0) CALL write_eig(chase_eig_id,ikpt,jsp,nev+nex,nev+nex,&
         eigenvalues(:(nev+nex)),zmat=zMatTemp)
381

382
    !Back-Transform
383
    IF (smat%l_real) THEN
384 385
       CALL PDTRTRI('U','N',smat%global_size1,smat%data_r,1,1,smat%blacs_desc,info)
       CALL PDGEMM('N','N',smat%global_size1,smat%global_size1,smat%global_size1,1.0,smat%data_r,1,1,smat%blacs_desc,zmatTemp%data_r,1,1,zmattemp%blacs_desc,0.0,hmat%data_r,1,1,hmat%blacs_desc)
386
    ELSE
387
       STOP 'chase no complex'
388
    ENDIF
389 390 391 392 393 394 395 396 397 398 399 400 401 402 403
!!$    CALL hmat%copy(zmatTemp,1,1) !Copy matrix into distributed form
!!$    call zmatTemp%free()
!!$    
!!$    ! --> recover the generalized eigenvectors z by solving z' = l^t * z
!!$    IF (smat%l_real) THEN
!!$       CALL pdtrtrs('U','N','N',hmat%global_size1,hmat%global_size1,smat%data_r,1,1,smat%blacs_desc,&
!!$            hmat%data_r,1,1,smat%blacs_desc,info)
!!$    ELSE
!!$       CALL pztrtrs('U','N','N',hmat%global_size1,hmat%global_size1,smat%data_c,1,1,smat%blacs_desc,&
!!$            hmat%data_c,1,1,smat%blacs_desc,info)
!!$    END IF
!!$    IF (info.NE.0) THEN
!!$       WRITE (6,*) 'Error in p?trtrs: info =',info
!!$       CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
!!$    ENDIF
404 405 406 407 408 409 410 411

    !     Redistribute eigvec from ScaLAPACK distribution to each process
    !     having all eigenvectors corresponding to his eigenvalues as above
    !
    ALLOCATE(t_mpimat::zmat)
    CALL zmat%init(hmat%l_real,hmat%global_size1,hmat%global_size1,hmat%mpi_com,.FALSE.)
    CALL zmat%copy(hmat,1,1)

Daniel Wortmann's avatar
Daniel Wortmann committed
412 413 414 415 416 417
    !Distribute eigenvalues over PEs
    ne=0
    DO i=myid+1,nev,np
       ne=ne+1
       eig(ne)=eigenvalues(i)
    ENDDO
418 419
  END SUBROUTINE chase_diag_MPI

420 421 422 423 424 425 426 427 428 429 430
  SUBROUTINE priv_init_chasempimat(hmat,mat,nev,nex)
    USE m_types_mpimat
    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
Daniel Wortmann's avatar
Daniel Wortmann committed
431
    INTEGER     :: k,i,j,l
432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472
    INTEGER     :: ierr

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

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

    mat%mpi_com=hmat%mpi_com
    mat%global_size1=hmat%global_size1
    mat%global_size2=hmat%global_size1
    mat%l_real=hmat%l_real
   
    !Determine rank and no of processors
    CALL MPI_COMM_RANK(hmat%mpi_com,myid,ierr)
    CALL MPI_COMM_SIZE(hmat%mpi_com,np,ierr)

    !Init ChASE
    IF (mat%l_real) THEN
       CALL mpi_dchase_init(hmat%mpi_com,mat%global_size1, nev, nex, rowoff,coloff,rowlen,collen,&
            mat%nprow,mat%npcol,myrow,mycol)
    ELSE
       CALL mpi_zchase_init(hmat%mpi_com,mat%global_size1, nev, nex, rowoff,coloff,rowlen,collen,&
            mat%nprow,mat%npcol,myrow,mycol)
    ENDIF

    !Determine block-sizes
    CALL MPI_ALLREDUCE(rowlen,nbr,1,MPI_INTEGER,MPI_MAX,mat%mpi_com,ierr)
    CALL MPI_ALLREDUCE(collen,nbc,1,MPI_INTEGER,MPI_MAX,mat%mpi_com,ierr)

    
    ALLOCATE(iusermap(mat%nprow,mat%npcol))
    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)
    CALL MPI_ALLREDUCE(MPI_IN_PLACE, iusermap, np,MPI_INTEGER,MPI_MAX,mat%mpi_com,ierr)
    !Get the Blacs default context
    CALL BLACS_GET(0,0,mat%blacs_ctext)
    ! Create the Grid
    CALL BLACS_GRIDMAP(mat%blacs_ctext,iusermap,mat%nprow,mat%nprow,mat%npcol)

Daniel Wortmann's avatar
Daniel Wortmann committed
473
   
474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499
    !Now create the matrix
    mat%matsize1=numroc(mat%global_size1,nbr,myrow,0,mat%nprow)
    mat%matsize2=numroc(mat%global_size1,nbc,mycol,0,mat%npcol)
    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
    CALL descinit(mat%blacs_desc,mat%global_size1,mat%global_size2,nbr,nbc,0,0,mat%blacs_ctext,mat%matsize1,ierr)
    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

500

501
#endif
502
  END MODULE m_chase_diag