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