LCOV - code coverage report
Current view: top level - diagonalization - chase_diag.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 5 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 2 0.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !
       6             : ! @authors: Miriam Hinzen, Gregor Michalicek
       7             : ! Added MPI implementation, DW 2018
       8             : !--------------------------------------------------------------------------------
       9             : MODULE m_chase_diag
      10             : #ifdef CPP_CHASE
      11             :   USE m_judft
      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,DIMENSION,input,atoms,kpts,noco,l_real)
      92             :     USE m_types_mpimat
      93             :     USE m_types_setup
      94             :     USE m_types_mpi
      95             :     USE m_judft
      96             :     USE m_eig66_io
      97             : 
      98             :     IMPLICIT NONE
      99             : 
     100             :     TYPE(t_mpi),               INTENT(IN)    :: mpi
     101             :     TYPE(t_dimension),         INTENT(IN)    :: dimension
     102             :     TYPE(t_input),             INTENT(IN)    :: input
     103             :     TYPE(t_atoms),             INTENT(IN)    :: atoms
     104             :     TYPE(t_kpts),              INTENT(IN)    :: kpts
     105             :     TYPE(t_noco),              INTENT(IN)    :: noco
     106             :     LOGICAL,                   INTENT(IN)    :: l_real
     107             : 
     108             :     INTEGER                                  :: nevd, nexd
     109             :     CHARACTER(len=1000)::arg
     110             :     
     111             :    
     112             :     scale_distance=1E-3
     113             :     !IF (judft_was_argument("-chase_tol_scale")) THEN
     114             :     !   arg=juDFT_string_for_argument("-chase_tol_scale")
     115             :     !   READ(arg,*) scale_distance
     116             :     !ENDIF
     117             :     
     118             :     IF (TRIM(juDFT_string_for_argument("-diag"))=="chase") THEN
     119             :        nevd = min(dimension%neigd,dimension%nvd+atoms%nlotot)
     120             :        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             :     END IF
     124             :   END SUBROUTINE init_chase
     125             : #endif
     126             :   
     127           0 :    SUBROUTINE chase_diag(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
     128             :     USE m_types_mpimat
     129             :     USE m_types_mat
     130             :     USE m_judft
     131             :     USE iso_c_binding
     132             :     USE m_eig66_io
     133             : 
     134             :     !Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
     135             :     IMPLICIT NONE
     136             : 
     137             :     CLASS(t_mat),              INTENT(INOUT) :: hmat,smat
     138             :     INTEGER,                   INTENT(IN)    :: ikpt
     139             :     INTEGER,                   INTENT(IN)    :: jsp
     140             :     INTEGER,                   INTENT(IN)    :: iter
     141             :     INTEGER,                   INTENT(INOUT) :: ne
     142             :     CLASS(t_mat), ALLOCATABLE, INTENT(OUT)   :: zmat
     143             :     REAL,                      INTENT(OUT)   :: eig(:)
     144             : #ifdef CPP_CHASE
     145             :     !Choose serial or parallel solver
     146             :     SELECT TYPE(hmat)
     147             :     CLASS is (t_mpimat)
     148             :        SELECT TYPE(smat)
     149             :        CLASS is (t_mpimat)
     150             :           CALL chase_diag_MPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
     151             :        CLASS default
     152             :           CALL judft_error("Inconsistent matrix setup")
     153             :        END SELECT
     154             :     CLASS is (t_mat)
     155             :        SELECT TYPE(smat)
     156             :        CLASS is (t_mat)
     157             :           CALL chase_diag_noMPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
     158             :        CLASS default
     159             :           CALL judft_error("Inconsistent matrix setup")
     160             :        END SELECT
     161             :     END SELECT
     162             : #endif    
     163           0 :   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             :     USE m_judft
     169             :     USE iso_c_binding
     170             :     USE m_eig66_io
     171             : 
     172             :     !Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
     173             :     IMPLICIT NONE
     174             : 
     175             :     TYPE(t_mat),               INTENT(INOUT) :: hmat,smat
     176             :     INTEGER,                   INTENT(IN)    :: ikpt
     177             :     INTEGER,                   INTENT(IN)    :: jsp
     178             :     INTEGER,                   INTENT(IN)    :: iter
     179             :     INTEGER,                   INTENT(INOUT) :: ne
     180             :     CLASS(t_mat), ALLOCATABLE, INTENT(OUT)   :: zmat
     181             :     REAL,                      INTENT(OUT)   :: eig(:)
     182             : 
     183             :     INTEGER            :: i, j, nev, nex, nbands
     184             :     INTEGER            :: info
     185             : 
     186             :     CLASS(t_Mat),              ALLOCATABLE  :: zMatTemp
     187             :     REAL(c_double),            ALLOCATABLE  :: eigenvalues(:)
     188             : 
     189             :     ALLOCATE(t_mat::zmat)
     190             :     CALL zmat%alloc(hmat%l_real,hmat%matsize1,ne)
     191             : 
     192             :     nev = min(ne,hmat%matsize1)
     193             :     nex = min(max(nev/4, 45), hmat%matsize1-nev) !dimensioning for workspace
     194             : 
     195             :     ALLOCATE(eigenvalues(nev+nex))
     196             :     eigenvalues = 0.0
     197             : 
     198             :     ALLOCATE(t_mat::zmatTemp)
     199             :     CALL zMatTemp%alloc(hmat%l_real,hmat%matsize1,nev+nex)
     200             : 
     201             :     IF (hmat%l_real) THEN
     202             : 
     203             :        ! --> start with Cholesky factorization of b ( so that b = l * l^t)
     204             :        ! --> b is overwritten by l
     205             :        CALL dpotrf('U',smat%matsize1,smat%data_r,SIZE(smat%data_r,1),info)
     206             :        IF (info.NE.0) THEN
     207             :           WRITE (*,*) 'Error in dpotrf: info =',info
     208             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     209             :        ENDIF
     210             : 
     211             :        ! --> now reduce a * z = eig * b * z to the standard form a' * z' = eig * z' 
     212             :        ! --> where a' = (l)^-1 * a * (l^t)^-1 and z' = l^t * z
     213             :        CALL dsygst(1,'U',smat%matsize1,hmat%data_r,SIZE(hmat%data_r,1),smat%data_r,SIZE(smat%data_r,1),info)
     214             :        IF (info.NE.0) THEN
     215             :           WRITE (6,*) 'Error in dsygst: info =',info
     216             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     217             :        ENDIF
     218             : 
     219             :        ! --> solve a' * z' = eig * z' for eigenvalues eig between lb und ub
     220             : 
     221             :        zMatTemp%data_r = 0.0
     222             : 
     223             :        do j = 1, hmat%matsize1
     224             :           do i = 1, j
     225             :              hmat%data_r(j,i) = hmat%data_r(i,j)
     226             :           end do
     227             :        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             :        ne = nev
     236             : 
     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             :        IF (info.NE.0) THEN
     243             :           WRITE (6,*) 'Error in dtrtrs: info =',info
     244             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     245             :        ENDIF
     246             : 
     247             :        DO i = 1, ne
     248             :           DO j = 1, hmat%matsize1
     249             :              zmat%data_r(j,i) = zMatTemp%data_r(j,i)
     250             :           END DO
     251             :           eig(i) = eigenvalues(i)
     252             :        END DO
     253             : 
     254             : 
     255             :     ELSE
     256             : 
     257             :        ! --> start with Cholesky factorization of b ( so that b = l * l^t)
     258             :        ! --> b is overwritten by l
     259             :        CALL zpotrf('U',smat%matsize1,smat%data_c,SIZE(smat%data_c,1),info)
     260             :        IF (info.NE.0) THEN
     261             :           WRITE (*,*) 'Error in zpotrf: info =',info
     262             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     263             :        ENDIF
     264             : 
     265             :        ! --> now reduce a * z = eig * b * z to the standard form a' * z' = eig * z' 
     266             :        ! --> where a' = (l)^-1 * a * (l^t)^-1 and z' = l^t * z
     267             :        CALL zhegst(1,'U',smat%matsize1,hmat%data_c,SIZE(hmat%data_c,1),smat%data_c,SIZE(smat%data_c,1),info)
     268             :        IF (info.NE.0) THEN
     269             :           WRITE (6,*) 'Error in zhegst: info =',info
     270             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     271             :        ENDIF
     272             : 
     273             :        ! --> solve a' * z' = eig * z' for eigenvalues eig between lb und ub
     274             : 
     275             :        zMatTemp%data_c = CMPLX(0.0,0.0)
     276             : 
     277             :        do j = 1, hmat%matsize1
     278             :           do i = 1, j
     279             :              hmat%data_c(j,i) = conjg(hmat%data_c(i,j))
     280             :           end do
     281             :        end do
     282             : 
     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             :        ne = nev
     291             : 
     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             :        IF (info.NE.0) THEN
     298             :           WRITE (6,*) 'Error in ztrtrs: info =',info
     299             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     300             :        ENDIF
     301             : 
     302             :        DO i = 1, ne
     303             :           DO j = 1, hmat%matsize1
     304             :              zmat%data_c(j,i) = zMatTemp%data_c(j,i)
     305             :           END DO
     306             :           eig(i) = eigenvalues(i)
     307             :        END DO
     308             : 
     309             :     ENDIF
     310             :     IF (info.NE.0) CALL judft_error("Diagonalization via ChASE failed", calledby = 'chase_diag')
     311             :   END SUBROUTINE chase_diag_noMPI
     312             : 
     313             :   SUBROUTINE chase_diag_MPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
     314             :     use m_types_mpimat
     315             :     USE m_types_mat
     316             :     USE m_judft
     317             :     USE iso_c_binding
     318             :     USE m_eig66_io
     319             :     
     320             :     !Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
     321             :     IMPLICIT NONE
     322             : 
     323             :     TYPE(t_mpimat),            INTENT(INOUT)    :: hmat,smat
     324             :     INTEGER,                   INTENT(IN)       :: ikpt
     325             :     INTEGER,                   INTENT(IN)       :: jsp
     326             :     INTEGER,                   INTENT(IN)       :: iter
     327             :     INTEGER,                   INTENT(INOUT)    :: ne
     328             :     CLASS(t_mat), ALLOCATABLE, INTENT(OUT)      :: zmat
     329             :     REAL,                      INTENT(OUT)      :: eig(:)
     330             : 
     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             :     TYPE(t_mpimat)                              :: chase_mat
     337             :     REAL,                          ALLOCATABLE  :: eigenvalues(:)
     338             : 
     339             :     REAL :: t1,t2,t3,t4
     340             :     include 'mpif.h'
     341             : 
     342             :     CALL CPU_TIME(t1)
     343             :     CALL MPI_COMM_RANK(hmat%blacsdata%mpi_com,myid,info)
     344             :     CALL MPI_COMM_SIZE(hmat%blacsdata%mpi_com,np,info)
     345             :     smat%blacsdata%blacs_desc=hmat%blacsdata%blacs_desc
     346             : 
     347             :     call smat%generate_full_matrix()
     348             :     call hmat%generate_full_matrix()
     349             :     !Transform to standard problem using SCALAPACK
     350             :     IF (hmat%l_real) THEN
     351             :        CALL pdpotrf('U',smat%global_size1,smat%data_r,1,1,smat%blacsdata%blacs_desc,info)
     352             :        IF (info.NE.0) THEN
     353             :           WRITE (*,*) 'Error in pdpotrf: info =',info
     354             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     355             :        ENDIF
     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             :         IF (info.NE.0) THEN
     359             :           WRITE (6,*) 'Error in pdsygst: info =',info
     360             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     361             :        ENDIF
     362             :     ELSE
     363             :        CALL pzpotrf('U',smat%global_size1,smat%data_c,1,1,smat%blacsdata%blacs_desc,info)
     364             :        IF (info.NE.0) THEN
     365             :           WRITE (*,*) 'Error in pzpotrf: info =',info
     366             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     367             :        ENDIF
     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             :         IF (info.NE.0) THEN
     371             :           WRITE (6,*) 'Error in pzhegst: info =',info
     372             :           CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     373             :        ENDIF
     374             :     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             :     ALLOCATE(eigenvalues(nev+nex))
     383             :     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             :     IF (hmat%l_real) THEN
     388             :        IF(iter.EQ.1) THEN
     389             :           CALL CPU_TIME(t2)
     390             :           CALL mpi_chase_r(chase_mat%data_r, zMatTemp%data_r, eigenvalues,  25, 1E-4, 'R', 'S' )
     391             :           CALL CPU_TIME(t3)
     392             :        ELSE
     393             :           CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
     394             :           CALL CPU_TIME(t2)
     395             :           CALL mpi_chase_r(chase_mat%data_r,  zMatTemp%data_r, eigenvalues, 25, tol, 'A', 'S' )
     396             :           CALL CPU_TIME(t3)
     397             :        END IF
     398             :     ELSE
     399             :        IF(iter.EQ.1) THEN
     400             :           CALL CPU_TIME(t2)
     401             :           CALL mpi_chase_c(chase_mat%data_c,  zMatTemp%data_c, eigenvalues,  25, 1E-4, 'R', 'S' )
     402             :           CALL CPU_TIME(t3)
     403             :        ELSE
     404             :           CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
     405             :           CALL CPU_TIME(t2)
     406             :           CALL mpi_chase_c(chase_mat%data_c,  zMatTemp%data_c, eigenvalues,  25, tol, 'A', 'S' )
     407             :           CALL CPU_TIME(t3)
     408             :        END IF
     409             :     ENDIF
     410             :     
     411             :     ne = nev
     412             :     IF (myid==0) CALL write_eig(chase_eig_id,ikpt,jsp,nev+nex,nev+nex,&
     413             :          eigenvalues(:(nev+nex)),zmat=zMatTemp)
     414             : 
     415             :     CALL hmat%from_non_dist(zmattemp)
     416             :     call zmatTemp%free()
     417             :     
     418             :     ! --> recover the generalized eigenvectors z by solving z' = l^t * z
     419             :     IF (smat%l_real) THEN
     420             :        CALL pdtrtrs('U','N','N',hmat%global_size1,hmat%global_size1,smat%data_r,1,1,smat%blacsdata%blacs_desc,&
     421             :             hmat%data_r,1,1,smat%blacsdata%blacs_desc,info)
     422             :     ELSE
     423             :        CALL pztrtrs('U','N','N',hmat%global_size1,hmat%global_size1,smat%data_c,1,1,smat%blacsdata%blacs_desc,&
     424             :             hmat%data_c,1,1,smat%blacsdata%blacs_desc,info)
     425             :     END IF
     426             :     IF (info.NE.0) THEN
     427             :        WRITE (6,*) 'Error in p?trtrs: info =',info
     428             :        CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
     429             :     ENDIF
     430             : 
     431             :     !     Redistribute eigvec from ScaLAPACK distribution to each process
     432             :     !     having all eigenvectors corresponding to his eigenvalues as above
     433             :     !
     434             :     ALLOCATE(t_mpimat::zmat)
     435             :     CALL zmat%init(hmat%l_real,hmat%global_size1,hmat%global_size1,hmat%blacsdata%mpi_com,.FALSE.)
     436             :     CALL zmat%copy(hmat,1,1)
     437             : 
     438             :     !Distribute eigenvalues over PEs
     439             :     ne=0
     440             :     DO i=myid+1,nev,np
     441             :        ne=ne+1
     442             :        eig(ne)=eigenvalues(i)
     443             :     ENDDO
     444             : 
     445             :     CALL CPU_TIME(t4)
     446             : 
     447             :     IF (myid==0) THEN
     448             :        PRINT *,"Chase Prep:",t2-t1
     449             :        PRINT *,"Chase Call:",t3-t2
     450             :        PRINT *,"Chase Post:",t4-t3
     451             :        PRINT *,"Chase Total:",t4-t1
     452             :     ENDIF
     453             : 
     454             :   END SUBROUTINE chase_diag_MPI
     455             : 
     456             :   SUBROUTINE priv_init_chasempimat(hmat,mat,nev,nex)
     457             :     USE m_types_mpimat
     458             :     USE m_types_mat
     459             :     IMPLICIT NONE
     460             :     TYPE(t_mpimat),INTENT(INOUT)::hmat,mat
     461             :     INTEGER,INTENT(IN)          :: nev,nex
     462             :     INTEGER::nbc,nbr
     463             : 
     464             :     INCLUDE 'mpif.h'
     465             :     INTEGER     :: myrow,mycol
     466             :     INTEGER     :: npblacs,np,myid
     467             :     INTEGER     :: rowlen,collen,rowoff,coloff
     468             :     INTEGER     :: k,i,j,l
     469             :     INTEGER     :: ierr
     470             : 
     471             :     INTEGER,ALLOCATABLE :: iblacsnums(:),ihelp(:),iusermap(:,:)
     472             : 
     473             :     EXTERNAL descinit, blacs_get
     474             :     EXTERNAL blacs_pinfo, blacs_gridinit
     475             :     INTEGER,EXTERNAL::numroc,indxl2g
     476             : 
     477             :     ALLOCATE(mat%blacsdata)
     478             :     mat%blacsdata%mpi_com=hmat%blacsdata%mpi_com
     479             :     mat%global_size1=hmat%global_size1
     480             :     mat%global_size2=hmat%global_size1
     481             :     mat%l_real=hmat%l_real
     482             :    
     483             :     !Determine rank and no of processors
     484             :     CALL MPI_COMM_RANK(hmat%blacsdata%mpi_com,myid,ierr)
     485             :     CALL MPI_COMM_SIZE(hmat%blacsdata%mpi_com,np,ierr)
     486             : 
     487             :     !Init ChASE
     488             :     IF (mat%l_real) THEN
     489             :        CALL mpi_dchase_init(hmat%blacsdata%mpi_com,mat%global_size1, nev, nex, rowoff,coloff,rowlen,collen,&
     490             :             mat%blacsdata%nprow,mat%blacsdata%npcol,myrow,mycol)
     491             :     ELSE
     492             :        CALL mpi_zchase_init(hmat%blacsdata%mpi_com,mat%global_size1, nev, nex, rowoff,coloff,rowlen,collen,&
     493             :             mat%blacsdata%nprow,mat%blacsdata%npcol,myrow,mycol)
     494             :     ENDIF
     495             : 
     496             :     !Determine block-sizes
     497             :     CALL MPI_ALLREDUCE(rowlen,nbr,1,MPI_INTEGER,MPI_MAX,mat%blacsdata%mpi_com,ierr)
     498             :     CALL MPI_ALLREDUCE(collen,nbc,1,MPI_INTEGER,MPI_MAX,mat%blacsdata%mpi_com,ierr)
     499             : 
     500             :     
     501             :     ALLOCATE(iusermap(mat%blacsdata%nprow,mat%blacsdata%npcol))
     502             :     iusermap=-2
     503             :     !Get BLACS ranks for all MPI ranks
     504             :     CALL BLACS_PINFO(iusermap(myrow+1,mycol+1),npblacs)  ! iamblacs = local process rank (e.g. myid)
     505             :     CALL MPI_ALLREDUCE(MPI_IN_PLACE, iusermap, np,MPI_INTEGER,MPI_MAX,mat%blacsdata%mpi_com,ierr)
     506             :     !Get the Blacs default context
     507             :     CALL BLACS_GET(0,0,mat%blacsdata%blacs_desc(2))
     508             :     ! Create the Grid
     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
     513             :     mat%matsize1=numroc(mat%global_size1,nbr,myrow,0,mat%blacsdata%nprow)
     514             :     mat%matsize2=numroc(mat%global_size1,nbc,mycol,0,mat%blacsdata%npcol)
     515             :     IF (mat%l_real) THEN
     516             :        ALLOCATE(mat%data_r(mat%matsize1,mat%matsize2))
     517             :     ELSE
     518             :        ALLOCATE(mat%data_c(mat%matsize1,mat%matsize2))
     519             :     END IF
     520             :     !Check for consistency
     521             :     IF (mat%matsize1.NE.rowlen.OR.mat%matsize2.NE.collen) THEN
     522             :        PRINT *,myid,"R:",mat%matsize1,rowlen,nbr
     523             :        PRINT *,myid,"C:",mat%matsize2,collen,nbc
     524             :        CALL judft_error("Distribution failed for chase")
     525             :     ENDIF
     526             : 
     527             :     !Create blacs descriptor for chase matrix
     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             :     IF (ierr /=0 ) CALL judft_error('Creation of BLACS descriptor failed')
     530             : 
     531             :     !Copy data from hmat
     532             :     CALL mat%copy(hmat,1,1)
     533             : 
     534             :     
     535             :     
     536             :   END SUBROUTINE priv_init_chasempimat
     537             : 
     538             : 
     539             : #endif
     540             :   END MODULE m_chase_diag

Generated by: LCOV version 1.13