Line data Source code
1 : ! (c) Copyright Pavel Kus, 2017, MPCDF
2 : !
3 : ! This file is part of ELPA.
4 : !
5 : ! The ELPA library was originally created by the ELPA consortium,
6 : ! consisting of the following organizations:
7 : !
8 : ! - Max Planck Computing and Data Facility (MPCDF), formerly known as
9 : ! Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
10 : ! - Bergische Universität Wuppertal, Lehrstuhl für angewandte
11 : ! Informatik,
12 : ! - Technische Universität München, Lehrstuhl für Informatik mit
13 : ! Schwerpunkt Wissenschaftliches Rechnen ,
14 : ! - Fritz-Haber-Institut, Berlin, Abt. Theorie,
15 : ! - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
16 : ! Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
17 : ! and
18 : ! - IBM Deutschland GmbH
19 : !
20 : !
21 : ! More information can be found here:
22 : ! http://elpa.mpcdf.mpg.de/
23 : !
24 : ! ELPA is free software: you can redistribute it and/or modify
25 : ! it under the terms of the version 3 of the license of the
26 : ! GNU Lesser General Public License as published by the Free
27 : ! Software Foundation.
28 : !
29 : ! ELPA is distributed in the hope that it will be useful,
30 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
31 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32 : ! GNU Lesser General Public License for more details.
33 : !
34 : ! You should have received a copy of the GNU Lesser General Public License
35 : ! along with ELPA. If not, see <http://www.gnu.org/licenses/>
36 : !
37 : ! ELPA reflects a substantial effort on the part of the original
38 : ! ELPA consortium, and we ask you to respect the spirit of the
39 : ! license that we chose: i.e., please contribute any changes you
40 : ! may have back to the original ELPA library distribution, and keep
41 : ! any derivatives of ELPA under the same license that we chose for
42 : ! the original distribution, the GNU Lesser General Public License.
43 :
44 : ! compute all eigenvectors
45 : subroutine solve_p&
46 : &BLAS_CHAR_AND_SY_OR_HE&
47 192 : &evd(na, a, sc_desc, ev, z)
48 : implicit none
49 : #include "../../src/general/precision_kinds.F90"
50 : integer(kind=ik), intent(in) :: na
51 : MATH_DATATYPE(kind=rck), intent(in) :: a(:,:)
52 : MATH_DATATYPE(kind=rck), intent(inout) :: z(:,:)
53 : real(kind=rk), intent(inout) :: ev(:)
54 : integer(kind=ik), intent(in) :: sc_desc(:)
55 : integer(kind=ik) :: info, lwork, liwork, lrwork
56 192 : MATH_DATATYPE(kind=rck), allocatable :: work(:)
57 192 : real(kind=rk), allocatable :: rwork(:)
58 192 : integer, allocatable :: iwork(:)
59 :
60 192 : allocate(work(1), iwork(1), rwork(1))
61 :
62 : ! query for required workspace
63 : #ifdef REALCASE
64 : call p&
65 : &BLAS_CHAR&
66 96 : &syevd('V', 'U', na, a, 1, 1, sc_desc, ev, z, 1, 1, sc_desc, work, -1, iwork, -1, info)
67 : #endif
68 : #ifdef COMPLEXCASE
69 : call p&
70 : &BLAS_CHAR&
71 96 : &heevd('V', 'U', na, a, 1, 1, sc_desc, ev, z, 1, 1, sc_desc, work, -1, rwork, -1, iwork, -1, info)
72 : #endif
73 : ! write(*,*) "computed sizes", lwork, liwork, "required sizes ", work(1), iwork(1)
74 192 : lwork = work(1)
75 192 : liwork = iwork(1)
76 192 : deallocate(work, iwork)
77 192 : allocate(work(lwork), stat = info)
78 192 : allocate(iwork(liwork), stat = info)
79 : #ifdef COMPLEXCASE
80 96 : lrwork = rwork(1)
81 96 : deallocate(rwork)
82 96 : allocate(rwork(lrwork), stat = info)
83 : #endif
84 : ! the actuall call to the method
85 : #ifdef REALCASE
86 : call p&
87 : &BLAS_CHAR&
88 96 : &syevd('V', 'U', na, a, 1, 1, sc_desc, ev, z, 1, 1, sc_desc, work, lwork, iwork, liwork, info)
89 : #endif
90 : #ifdef COMPLEXCASE
91 : call p&
92 : &BLAS_CHAR&
93 96 : &heevd('V', 'U', na, a, 1, 1, sc_desc, ev, z, 1, 1, sc_desc, work, lwork, rwork, lrwork, iwork, liwork, info)
94 : #endif
95 :
96 192 : deallocate(iwork, work, rwork)
97 192 : end subroutine
98 :
99 :
100 : ! compute part of eigenvectors
101 : subroutine solve_p&
102 : &BLAS_CHAR_AND_SY_OR_HE&
103 192 : &evr(na, a, sc_desc, nev, ev, z)
104 : implicit none
105 : #include "../../src/general/precision_kinds.F90"
106 : integer(kind=ik), intent(in) :: na, nev
107 : MATH_DATATYPE(kind=rck), intent(in) :: a(:,:)
108 : MATH_DATATYPE(kind=rck), intent(inout) :: z(:,:)
109 : real(kind=rk), intent(inout) :: ev(:)
110 : integer(kind=ik), intent(in) :: sc_desc(:)
111 : integer(kind=ik) :: info, lwork, liwork, lrwork
112 192 : MATH_DATATYPE(kind=rck), allocatable :: work(:)
113 192 : real(kind=rk), allocatable :: rwork(:)
114 192 : integer, allocatable :: iwork(:)
115 : integer(kind=ik) :: comp_eigenval, comp_eigenvec, smallest_ev_idx, largest_ev_idx
116 :
117 192 : allocate(work(1), iwork(1), rwork(1))
118 192 : smallest_ev_idx = 1
119 192 : largest_ev_idx = nev
120 : ! query for required workspace
121 : #ifdef REALCASE
122 : call p&
123 : &BLAS_CHAR&
124 : &syevr('V', 'I', 'U', na, a, 1, 1, sc_desc, 0.0_rk, 0.0_rk, smallest_ev_idx, largest_ev_idx, &
125 96 : comp_eigenval, comp_eigenvec, ev, z, 1, 1, sc_desc, work, -1, iwork, -1, info)
126 : #endif
127 : #ifdef COMPLEXCASE
128 : call p&
129 : &BLAS_CHAR&
130 : &heevr('V', 'I', 'U', na, a, 1, 1, sc_desc, 0.0_rk, 0.0_rk, smallest_ev_idx, largest_ev_idx, &
131 96 : comp_eigenval, comp_eigenvec, ev, z, 1, 1, sc_desc, work, -1, rwork, -1, iwork, -1, info)
132 : #endif
133 : ! write(*,*) "computed sizes", lwork, liwork, "required sizes ", work(1), iwork(1)
134 192 : lwork = work(1)
135 192 : liwork = iwork(1)
136 192 : deallocate(work, iwork)
137 192 : allocate(work(lwork), stat = info)
138 192 : allocate(iwork(liwork), stat = info)
139 : #ifdef COMPLEXCASE
140 96 : lrwork = rwork(1)
141 96 : deallocate(rwork)
142 96 : allocate(rwork(lrwork), stat = info)
143 : #endif
144 : ! the actuall call to the method
145 : #ifdef REALCASE
146 : call p&
147 : &BLAS_CHAR&
148 : &syevr('V', 'I', 'U', na, a, 1, 1, sc_desc, 0.0_rk, 0.0_rk, smallest_ev_idx, largest_ev_idx, &
149 96 : comp_eigenval, comp_eigenvec, ev, z, 1, 1, sc_desc, work, lwork, iwork, liwork, info)
150 : #endif
151 : #ifdef COMPLEXCASE
152 : call p&
153 : &BLAS_CHAR&
154 : &heevr('V', 'I', 'U', na, a, 1, 1, sc_desc, 0.0_rk, 0.0_rk, smallest_ev_idx, largest_ev_idx, &
155 96 : comp_eigenval, comp_eigenvec, ev, z, 1, 1, sc_desc, work, lwork, rwork, lrwork, iwork, liwork, info)
156 : #endif
157 192 : assert(comp_eigenval == nev)
158 192 : assert(comp_eigenvec == nev)
159 192 : deallocate(iwork, work, rwork)
160 192 : end subroutine
161 :
|
