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_lapack_singlePrec_diag
9 : USE m_judft
10 : USE m_constants
11 : IMPLICIT NONE
12 : PRIVATE
13 : integer,parameter:: sp=selected_real_kind(6)
14 :
15 : PUBLIC lapack_singlePrec_diag
16 :
17 : CONTAINS
18 :
19 0 : SUBROUTINE lapack_singlePrec_diag(hmat,smat,ne,eig,zmat)
20 :
21 : USE m_types_mat
22 : USE m_judft
23 :
24 : !Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
25 : IMPLICIT NONE
26 :
27 : TYPE(t_mat), INTENT(INOUT) :: hmat,smat
28 : INTEGER, INTENT(INOUT) :: ne
29 : CLASS(t_mat), ALLOCATABLE, INTENT(OUT) :: zmat
30 : REAL, INTENT(OUT) :: eig(:)
31 :
32 : INTEGER :: nev,lwork,liwork
33 : INTEGER :: info
34 :
35 :
36 :
37 0 : ALLOCATE(t_mat::zmat)
38 0 : CALL zmat%alloc(hmat%l_real,hmat%matsize1,ne)
39 :
40 :
41 0 : IF (hmat%l_real) THEN
42 :
43 : ! --> start with Cholesky factorization of b ( so that b = l * l^t)
44 : ! --> b is overwritten by l
45 0 : CALL dpotrf('U',smat%matsize1,smat%data_r,SIZE(smat%data_r,1),info)
46 0 : IF (info.NE.0) THEN
47 0 : WRITE (*,*) 'Error in dpotrf: info =',info
48 0 : CALL juDFT_error("Diagonalization failed",calledby="lapack_singlePrec_diag")
49 : ENDIF
50 :
51 : ! --> now reduce a * z = eig * b * z to the standard form a' * z' = eig * z'
52 : ! --> where a' = (l)^-1 * a * (l^t)^-1 and z' = l^t * z
53 0 : CALL dsygst(1,'U',smat%matsize1,hmat%data_r,SIZE(hmat%data_r,1),smat%data_r,SIZE(smat%data_r,1),info)
54 0 : IF (info.NE.0) THEN
55 0 : WRITE (oUnit,*) 'Error in dsygst: info =',info
56 0 : CALL juDFT_error("Diagonalization failed",calledby="lapack_singlePrec_diag")
57 : ENDIF
58 :
59 : ! --> solve a' * z' = eig * z' for eigenvalues eig between lb und ub
60 0 : BLOCK
61 0 : REAL(kind=sp),allocatable:: h(:,:),z(:,:),eigval(:),work(:)
62 0 : integer,allocatable :: iwork(:),ifail(:)
63 0 : Allocate(h(size(hmat%data_r,1),size(hmat%data_r,2)))
64 0 : Allocate(eigval(size(hmat%data_r,1)),ifail(size(hmat%data_r,1)))
65 0 : Allocate(z(size(hmat%data_r,1),ne))
66 0 : h=hmat%data_r
67 :
68 0 : allocate(work(1),iwork(1))
69 0 : call ssyevx('V','I','U',size(h,1),h,size(h,1),0.0,0.0,1,ne,0.0,nev,eigval,z,size(z,1),work,-1,iwork,liwork,ifail,info)
70 0 : lwork=work(1)
71 0 : liwork=iwork(1)
72 0 : deallocate(work,iwork)
73 0 : allocate(work(lwork),iwork(liwork))
74 :
75 0 : call ssyevx('V','I','U',size(h,1),h,size(h,1),0.0,0.0,1,ne,0.0,nev,eigval,z,size(z,1),work,lwork,iwork,liwork,ifail,info)
76 :
77 0 : eig=eigval(:ne)
78 0 : zmat%data_r=z(:,:ne)
79 0 : deallocate(h,z,eigval,work,iwork)
80 : END BLOCK
81 : ! --> recover the generalized eigenvectors z by solving z' = l^t * z
82 0 : CALL dtrtrs('U','N','N',hmat%matsize1,nev,smat%data_r,smat%matsize1,zMat%data_r,zmat%matsize1,info)
83 0 : IF (info.NE.0) THEN
84 0 : WRITE (oUnit,*) 'Error in dtrtrs: info =',info
85 0 : CALL juDFT_error("Diagonalization failed",calledby="lapack_singlePrec_diag")
86 : ENDIF
87 :
88 :
89 : ELSE
90 :
91 : ! --> start with Cholesky factorization of b ( so that b = l * l^t)
92 : ! --> b is overwritten by l
93 0 : CALL zpotrf('U',smat%matsize1,smat%data_c,SIZE(smat%data_c,1),info)
94 0 : IF (info.NE.0) THEN
95 0 : WRITE (*,*) 'Error in zpotrf: info =',info
96 0 : CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
97 : ENDIF
98 :
99 : ! --> now reduce a * z = eig * b * z to the standard form a' * z' = eig * z'
100 : ! --> where a' = (l)^-1 * a * (l^t)^-1 and z' = l^t * z
101 0 : CALL zhegst(1,'U',smat%matsize1,hmat%data_c,SIZE(hmat%data_c,1),smat%data_c,SIZE(smat%data_c,1),info)
102 0 : IF (info.NE.0) THEN
103 0 : WRITE (oUnit,*) 'Error in zhegst: info =',info
104 0 : CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
105 : ENDIF
106 :
107 : ! --> solve a' * z' = eig * z' for eigenvalues eig between lb und ub
108 0 : BLOCK
109 0 : COMPLEX(kind=sp),allocatable:: h(:,:),z(:,:),work(:)
110 0 : REAL(kind=sp),allocatable:: eigval(:),rwork(:)
111 0 : integer,allocatable :: iwork(:),ifail(:)
112 0 : Allocate(h(size(hmat%data_c,1),size(hmat%data_c,2)))
113 0 : Allocate(eigval(size(hmat%data_c,1)),ifail(size(hmat%data_c,1)))
114 0 : Allocate(z(size(hmat%data_c,1),ne),rwork(7*size(hmat%data_c,1)))
115 0 : h=hmat%data_c
116 :
117 0 : allocate(work(1),iwork(5*size(hmat%data_c,1)))
118 0 : call cheevx('V','I','U',size(h,1),h,size(h,1),0.0,0.0,1,ne,0.0,nev,eigval,z,size(z,1),work,-1,rwork,iwork,ifail,info)
119 0 : lwork=work(1)
120 0 : deallocate(work)
121 0 : allocate(work(lwork))
122 :
123 0 : call cheevx('V','I','U',size(h,1),h,size(h,1),0.0,0.0,1,ne,0.0,nev,eigval,z,size(z,1),work,lwork,rwork,iwork,ifail,info)
124 0 : eig=eigval(:ne)
125 0 : zmat%data_c=z(:,:ne)
126 0 : deallocate(h,z,eigval,work,rwork,iwork)
127 : END BLOCK
128 :
129 : ! --> recover the generalized eigenvectors z by solving z' = l^t * z
130 0 : CALL ztrtrs('U','N','N',hmat%matsize1,nev,smat%data_c,smat%matsize1,zMat%data_c,zmat%matsize1,info)
131 0 : IF (info.NE.0) THEN
132 0 : WRITE (oUnit,*) 'Error in ztrtrs: info =',info
133 0 : CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
134 : ENDIF
135 :
136 :
137 : ENDIF
138 0 : IF (info.NE.0) CALL judft_error("Diagonalization via ChASE failed", calledby = 'chase_diag')
139 0 : END SUBROUTINE lapack_singlePREC_diag
140 :
141 : END MODULE m_lapack_singlePrec_diag
|
