1!--------------------------------------------------------------------------------------------------!
2!  DFTB+: general package for performing fast atomistic simulations                                !
3!  Copyright (C) 2006 - 2019  DFTB+ developers group                                               !
4!                                                                                                  !
5!  See the LICENSE file for terms of usage and distribution.                                       !
6!--------------------------------------------------------------------------------------------------!
7
8#:include 'common.fypp'
9
10!> Contains F90 wrapper functions for some commonly used lapack calls needed in the code.
11!> Contains some fixes for lapack 3.0 bugs, if this gets corrected in lapack 4.x they should be
12!> removed.
13module dftbp_eigensolver
14  use dftbp_assert
15  use dftbp_message
16  use dftbp_accuracy, only : rsp, rdp
17  use dftbp_blas
18  use dftbp_lapack
19#:if WITH_GPU
20  use magma
21#:endif
22  implicit none
23  private
24
25  public :: heev, hegv, hegvd, gvr, bgv
26#:if WITH_GPU
27  public :: gpu_gvd
28#:endif
29
30
31  !> Used to return runtime diagnostics
32  character(len=100) :: error_string
33
34
35  !> Simple eigensolver for a symmetric/Hermitian matrix
36  !> Caveat: the matrix a is overwritten
37  interface heev
38    module procedure real_ssyev
39    module procedure dble_dsyev
40    module procedure cmplx_cheev
41    module procedure dblecmplx_zheev
42  end interface heev
43
44
45  !> Simple eigensolver for a symmetric/Hermitian generalized matrix problem
46  !> caveat: the matrix a is overwritten
47  !> caveat: the matrix b is overwritten with Cholesky factorization
48  interface hegv
49    module procedure real_ssygv
50    module procedure dble_dsygv
51    module procedure cmplx_chegv
52    module procedure dblecmplx_zhegv
53  end interface hegv
54
55
56  !> Simple eigensolver for a symmetric/Hermitian generalized matrix problem using divide and
57  !> conquer eigensolver
58  !> caveat: the matrix a is overwritten
59  !> caveat: the matrix b is overwritten with Cholesky factorization
60  interface hegvd
61    module procedure real_ssygvd
62    module procedure dble_dsygvd
63    module procedure cmplx_chegvd
64    module procedure dblecmplx_zhegvd
65  end interface hegvd
66
67
68  !> Simple eigensolver for a symmetric/Hermitian generalized matrix problem using the lapack
69  !> relatively robust representation solver, based on the SYGV source. If the requested number of
70  !> eigenvalues is lower than the size of H/S suspace mode is used (optionally the range can be set
71  !> using il and ul) to return the lowest eigenvalues/vectors of number size(w)
72  interface gvr
73    module procedure real_ssygvr
74    module procedure dble_dsygvr
75    module procedure cmplx_chegvr
76    module procedure dblecmplx_zhegvr
77  end interface
78
79
80  !> Eigensolver for a symmetric/Hermitian banded generalized matrix
81  !> problem of the form A*x=(lambda)*B*x
82  interface bgv
83    module procedure real_ssbgv
84    module procedure dble_dsbgv
85    module procedure cmplx_chbgv
86    module procedure dblecmplx_zhbgv
87  end interface
88
89#:if WITH_GPU
90  !> Divide and conquer MAGMA GPU eigensolver
91  interface gpu_gvd
92    module procedure real_magma_ssygvd
93    module procedure dble_magma_dsygvd
94    module procedure cmplx_magma_chegvd
95    module procedure dblecmplx_magma_zhegvd
96  end interface
97#:endif
98
99contains
100
101  !> Real eigensolver for a symmetric matrix
102  subroutine real_ssyev(a,w,uplo,jobz)
103
104    !> contains the matrix for the solver, returns as eigenvectors if requested (matrix always
105    !> overwritten on return anyway)
106    real(rsp), intent(inout) :: a(:,:)
107
108    !> eigenvalues
109    real(rsp), intent(out) :: w(:)
110
111    !> upper or lower triangle of the matrix
112    character, intent(in) :: uplo
113
114    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
115    character, intent(in) :: jobz
116
117    real(rsp), allocatable :: work(:)
118    integer n, info
119    integer :: int_idealwork
120    real(rsp) :: idealwork(1)
121    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
122    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
123    @:ASSERT(all(shape(a)==size(w,dim=1)))
124    n=size(a,dim=1)
125    @:ASSERT(n>0)
126    call ssyev(jobz, uplo, n, a, n, w, idealwork, -1, info)
127    if (info/=0) then
128      call error("Failure in SSYEV to determine optimum workspace")
129    endif
130    int_idealwork=floor(idealwork(1))
131    allocate(work(int_idealwork))
132    call ssyev(jobz, uplo, n, a, n, w, work, int_idealwork, info)
133    if (info/=0) then
134      if (info<0) then
13599000 format ('Failure in diagonalisation routine ssyev,', &
136          & ' illegal argument at position ',i6)
137        write(error_string, 99000) info
138        call error(error_string)
139      else
14099010 format ('Failure in diagonalisation routine ssyev,', &
141          & ' diagonal element ',i6,' did not converge to zero.')
142        write(error_string, 99010) info
143        call error(error_string)
144      endif
145    endif
146
147  end subroutine real_ssyev
148
149
150  !> Double precision eigensolver for a symmetric matrix
151  subroutine dble_dsyev(a,w,uplo,jobz)
152
153    !> contains the matrix for the solver, returns as eigenvectors if requested (matrix always
154    !> overwritten on return anyway)
155    real(rdp), intent(inout) :: a(:,:)
156
157    !> eigenvalues
158    real(rdp), intent(out) :: w(:)
159
160    !> upper or lower triangle of the matrix
161    character, intent(in) :: uplo
162
163    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
164    character, intent(in) :: jobz
165
166    real(rdp), allocatable :: work(:)
167    integer n, info
168    integer :: int_idealwork
169    real(rdp) :: idealwork(1)
170    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
171    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
172    @:ASSERT(all(shape(a)==size(w,dim=1)))
173    n=size(a,dim=1)
174    @:ASSERT(n>0)
175    call dsyev(jobz, uplo, n, a, n, w, idealwork, -1, info)
176    if (info/=0) then
177      call error("Failure in DSYEV to determine optimum workspace")
178    endif
179    int_idealwork=floor(idealwork(1))
180    allocate(work(int_idealwork))
181    call dsyev(jobz, uplo, n, a, n, w, work, int_idealwork, info)
182    if (info/=0) then
183      if (info<0) then
18499020 format ('Failure in diagonalisation routine dsyev,', &
185          & ' illegal argument at position ',i6)
186        write(error_string, 99020) info
187        call error(error_string)
188      else
18999030 format ('Failure in diagonalisation routine dsyev,', &
190          & ' diagonal element ',i6,' did not converge to zero.')
191        write(error_string, 99030) info
192        call error(error_string)
193      endif
194    endif
195
196  end subroutine dble_dsyev
197
198
199  !> Complex eigensolver for a Hermitian matrix
200  subroutine cmplx_cheev(a,w,uplo,jobz)
201
202    !> contains the matrix for the solver, returns as eigenvectors if requested (matrix always
203    !> overwritten on return anyway)
204    complex(rsp), intent(inout) :: a(:,:)
205
206    !> eigenvalues
207    real(rsp), intent(out) :: w(:)
208
209    !> upper or lower triangle of the matrix
210    character, intent(in) :: uplo
211
212    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
213    character, intent(in) :: jobz
214
215    real(rsp), allocatable :: rwork(:)
216    complex(rsp), allocatable :: work(:)
217    integer n, info
218    integer :: int_idealwork
219    complex(rsp) :: idealwork(1)
220    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
221    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
222    @:ASSERT(all(shape(a)==size(w,dim=1)))
223    n=size(a,dim=1)
224    @:ASSERT(n>0)
225    allocate(rwork(3*n-2))
226    call cheev(jobz, uplo, n, a, n, w, idealwork, -1, rwork, info)
227    if (info/=0) then
228      call error("Failure in CHEEV to determine optimum workspace")
229    endif
230    int_idealwork=floor(real(idealwork(1)))
231    allocate(work(int_idealwork))
232    call cheev(jobz, uplo, n, a, n, w, work, int_idealwork, rwork, info)
233    if (info/=0) then
234      if (info<0) then
23599040 format ('Failure in diagonalisation routine cheev,', &
236          & ' illegal argument at position ',i6)
237        write(error_string, 99040) info
238        call error(error_string)
239      else
24099050 format ('Failure in diagonalisation routine cheev,', &
241          & ' diagonal element ',i6,' did not converge to zero.')
242        write(error_string, 99050) info
243        call error(error_string)
244      endif
245    endif
246
247  end subroutine cmplx_cheev
248
249
250  !> Double complex eigensolver for a Hermitian matrix
251  subroutine dblecmplx_zheev(a,w,uplo,jobz)
252
253    !> contains the matrix for the solver, returns as eigenvectors if requested (matrix always
254    !> overwritten on return anyway)
255    complex(rdp), intent(inout) :: a(:,:)
256
257    !> eigenvalues
258    real(rdp), intent(out) :: w(:)
259
260    !> upper or lower triangle of the matrix
261    character, intent(in) :: uplo
262
263    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
264    character, intent(in) :: jobz
265
266    real(rdp), allocatable :: rwork(:)
267    complex(rdp), allocatable :: work(:)
268    integer n, info
269    integer :: int_idealwork
270    complex(rdp) :: idealwork(1)
271    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
272    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
273    @:ASSERT(all(shape(a)==size(w,dim=1)))
274    n=size(a,dim=1)
275    @:ASSERT(n>0)
276    allocate(rwork(3*n-2))
277    call zheev(jobz, uplo, n, a, n, w, idealwork, -1, rwork, info)
278    if (info/=0) then
279      call error("Failure in ZHEEV to determine optimum workspace")
280    endif
281    int_idealwork=floor(real(idealwork(1)))
282    allocate(work(int_idealwork))
283    call zheev(jobz, uplo, n, a, n, w, work, int_idealwork, rwork, info)
284    if (info/=0) then
285      if (info<0) then
28699060 format ('Failure in diagonalisation routine zheev,', &
287          & ' illegal argument at position ',i6)
288        write(error_string, 99060) info
289        call error(error_string)
290      else
29199070 format ('Failure in diagonalisation routine zheev,', &
292          & ' diagonal element ',i6,' did not converge to zero.')
293        write(error_string, 99070) info
294        call error(error_string)
295      endif
296    endif
297
298  end subroutine dblecmplx_zheev
299
300
301  !> Real eigensolver for generalized symmetric matrix problem
302  subroutine real_ssygv(a,b,w,uplo,jobz,itype)
303
304    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
305    !> overwritten on return anyway)
306    real(rsp), intent(inout) :: a(:,:)
307
308    !> contains the second matrix for the solver (overwritten by Cholesky factorization)
309    real(rsp), intent(inout) :: b(:,:)
310
311    !> eigenvalues
312    real(rsp), intent(out) :: w(:)
313
314    !> upper or lower triangle of both matrices
315    character, intent(in) :: uplo
316
317    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
318    character, intent(in) :: jobz
319
320    !> specifies the problem type to be solved 1:A*x=(lambda)*B*x, 2:A*B*x=(lambda)*x,
321    !> 3:B*A*x=(lambda)*x default is 1
322    integer, optional, intent(in) :: itype
323
324    real(rsp), allocatable :: work(:)
325    integer n, info, iitype
326    integer :: int_idealwork
327    real(rsp) :: idealwork(1)
328    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
329    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
330    @:ASSERT(all(shape(a)==shape(b)))
331    @:ASSERT(all(shape(a)==size(w,dim=1)))
332    n=size(a,dim=1)
333    @:ASSERT(n>0)
334    if (present(itype)) then
335      iitype = itype
336    else
337      iitype = 1
338    end if
339    @:ASSERT(iitype >= 1 .and. iitype <= 3 )
340    call ssygv(iitype, jobz, uplo, n, a, n, b, n, w, idealwork, -1, info)
341    if (info/=0) then
342       call error("Failure in SSYGV to determine optimum workspace")
343    endif
344    int_idealwork=floor(idealwork(1))
345    allocate(work(int_idealwork))
346    call ssygv(iitype, jobz, uplo, n, a, n, b, n, w, work, int_idealwork, info)
347    if (info/=0) then
348       if (info<0) then
34999160 format ('Failure in diagonalisation routine ssygv,', &
350          & ' illegal argument at position ',i6)
351          write(error_string, 99160) info
352          call error(error_string)
353       else if (info <= n) then
35499170 format ('Failure in diagonalisation routine ssygv,', &
355          & ' diagonal element ',i6,' did not converge to zero.')
356          write(error_string, 99170) info
357          call error(error_string)
358       else
35999180 format ('Failure in diagonalisation routine ssygv,', &
360          & ' non-positive definite overlap! Minor ',i6,' responsible.')
361          write(error_string, 99180) info - n
362          call error(error_string)
363       endif
364    endif
365
366  end subroutine real_ssygv
367
368
369  !> Double precision eigensolver for generalized symmetric matrix problem
370  subroutine dble_dsygv(a,b,w,uplo,jobz,itype)
371
372    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
373    !> overwritten on return anyway)
374    real(rdp), intent(inout) :: a(:,:)
375
376    !> contains the second matrix for the solver (overwritten by Cholesky factorization)
377    real(rdp), intent(inout) :: b(:,:)
378
379    !> eigenvalues
380    real(rdp), intent(out) :: w(:)
381
382    !> upper or lower triangle of both matrices
383    character, intent(in) :: uplo
384
385    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
386    character, intent(in) :: jobz
387
388    !> specifies the problem type to be solved 1:A*x=(lambda)*B*x, 2:A*B*x=(lambda)*x,
389    !> 3:B*A*x=(lambda)*x default is 1
390    integer, optional, intent(in) :: itype
391
392    real(rdp), allocatable :: work(:)
393    integer n, info, iitype
394    integer :: int_idealwork
395    real(rdp) :: idealwork(1)
396    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
397    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
398    @:ASSERT(all(shape(a)==shape(b)))
399    @:ASSERT(all(shape(a)==size(w,dim=1)))
400    n=size(a,dim=1)
401    @:ASSERT(n>0)
402    if (present(itype)) then
403      iitype = itype
404    else
405      iitype = 1
406    end if
407    @:ASSERT(iitype >= 1 .and. iitype <= 3 )
408    call dsygv(iitype, jobz, uplo, n, a, n, b, n, w, idealwork, -1, info)
409    if (info/=0) then
410       call error("Failure in DSYGV to determine optimum workspace")
411    endif
412    int_idealwork=floor(idealwork(1))
413    allocate(work(int_idealwork))
414    call dsygv(iitype, jobz, uplo, n, a, n, b, n, w, work, int_idealwork, info)
415    if (info/=0) then
416       if (info<0) then
41799190 format ('Failure in diagonalisation routine dsygv,', &
418          & ' illegal argument at position ',i6)
419          write(error_string, 99190) info
420          call error(error_string)
421       else if (info <= n) then
42299200 format ('Failure in diagonalisation routine dsygv,', &
423          & ' diagonal element ',i6,' did not converge to zero.')
424          write(error_string, 99200) info
425          call error(error_string)
426       else
42799210 format ('Failure in diagonalisation routine dsygv,', &
428          & ' non-positive definite overlap! Minor ',i6,' responsible.')
429          write(error_string, 99210) info - n
430          call error(error_string)
431       endif
432    endif
433
434  end subroutine dble_dsygv
435
436
437  !> Complex eigensolver for generalized Hermitian matrix problem
438  subroutine cmplx_chegv(a,b,w,uplo,jobz,itype)
439
440    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
441    !> overwritten on return anyway)
442    complex(rsp), intent(inout) :: a(:,:)
443
444    !> contains the second matrix for the solver (overwritten by Cholesky factorization)
445    complex(rsp), intent(inout) :: b(:,:)
446
447    !> eigenvalues
448    real(rsp), intent(out) :: w(:)
449
450    !> upper or lower triangle of both matrices
451    character, intent(in) :: uplo
452
453    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
454    character, intent(in) :: jobz
455
456    !> specifies the problem type to be solved 1:A*x=(lambda)*B*x, 2:A*B*x=(lambda)*x,
457    !> 3:B*A*x=(lambda)*x default is 1
458    integer, optional, intent(in) :: itype
459
460    complex(rsp), allocatable :: work(:)
461    real(rsp), allocatable :: rwork(:)
462    integer n, info, iitype, int_idealwork
463    complex(rsp) :: idealwork(1)
464
465    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
466    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
467    @:ASSERT(all(shape(a)==shape(b)))
468    @:ASSERT(all(shape(a)==size(w,dim=1)))
469    n=size(a,dim=1)
470    @:ASSERT(n>0)
471    if (present(itype)) then
472      iitype = itype
473    else
474      iitype = 1
475    end if
476    @:ASSERT(iitype >= 1 .and. iitype <= 3 )
477    allocate(rwork(3*n-2))
478    call CHEGV(iitype, jobz, uplo, n, a, n, b, n, w, idealwork, -1, rwork, info)
479    if (info/=0) then
480       call error("Failure in CHEGV to determine optimum workspace")
481    endif
482    int_idealwork=floor(real(idealwork(1)))
483    allocate(work(int_idealwork))
484    ! A*x = (lambda)*B*x upper triangles to be used
485    call chegv(iitype, 'V', 'L', n, a, n, b, n, w, work, int_idealwork, rwork, info)
486    if (info/=0) then
487       if (info<0) then
48899220 format ('Failure in diagonalisation routine chegv,', &
489          & ' illegal argument at position ',i6)
490          write(error_string, 99220) info
491          call error(error_string)
492       else if (info <= n) then
49399230 format ('Failure in diagonalisation routine chegv,', &
494          & ' diagonal element ',i6,' did not converge to zero.')
495          write(error_string, 99230) info
496          call error(error_string)
497       else
49899240 format ('Failure in diagonalisation routine chegv,', &
499          & ' non-positive definite overlap! Minor ',i6,' responsible.')
500          write(error_string, 99240) info - n
501          call error(error_string)
502       endif
503    endif
504
505  end subroutine cmplx_chegv
506
507
508  !> Double complex eigensolver for generalized Hermitian matrix problem
509  subroutine dblecmplx_zhegv(a,b,w,uplo,jobz,itype)
510
511    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
512    !> overwritten on return anyway)
513    complex(rdp), intent(inout) :: a(:,:)
514
515    !> contains the second matrix for the solver (overwritten by Cholesky factorization)
516    complex(rdp), intent(inout) :: b(:,:)
517
518    !> eigenvalues
519    real(rdp), intent(out) :: w(:)
520
521    !> upper or lower triangle of both matrices
522    character, intent(in) :: uplo
523
524    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
525    character, intent(in) :: jobz
526
527    !> specifies the problem type to be solved 1:A*x=(lambda)*B*x, 2:A*B*x=(lambda)*x,
528    !> 3:B*A*x=(lambda)*x default is 1
529    integer, optional, intent(in) :: itype
530
531    complex(rdp), allocatable :: work(:)
532    real(rdp), allocatable :: rwork(:)
533    integer n, info, iitype, int_idealwork
534    complex(rdp) :: idealwork(1)
535
536    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
537    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
538    @:ASSERT(all(shape(a)==shape(b)))
539    @:ASSERT(all(shape(a)==size(w,dim=1)))
540    n=size(a,dim=1)
541    @:ASSERT(n>0)
542    if (present(itype)) then
543      iitype = itype
544    else
545      iitype = 1
546    end if
547    @:ASSERT(iitype >= 1 .and. iitype <= 3 )
548    allocate(rwork(3*n-2))
549    call ZHEGV(iitype, jobz, uplo, n, a, n, b, n, w, idealwork, -1, rwork, info)
550    if (info/=0) then
551       call error("Failure in CHEGV to determine optimum workspace")
552    endif
553    int_idealwork=floor(real(idealwork(1)))
554    allocate(work(int_idealwork))
555    ! A*x = (lambda)*B*x upper triangles to be used
556    call zhegv(iitype, 'V', 'L', n, a, n, b, n, w, work, int_idealwork, rwork, info)
557    if (info/=0) then
558       if (info<0) then
55999250 format ('Failure in diagonalisation routine zhegv,', &
560          & ' illegal argument at position ',i6)
561          write(error_string, 99250) info
562          call error(error_string)
563       else if (info <= n) then
56499260 format ('Failure in diagonalisation routine zhegv,', &
565          & ' diagonal element ',i6,' did not converge to zero.')
566          write(error_string, 99260) info
567          call error(error_string)
568       else
56999270 format ('Failure in diagonalisation routine zhegv,', &
570          & ' non-positive definite overlap! Minor ',i6,' responsible.')
571          write(error_string, 99270) info - n
572          call error(error_string)
573       endif
574    endif
575
576  end subroutine dblecmplx_zhegv
577
578
579  !> Real eigensolver for generalized symmetric matrix problem - divide and conquer
580  subroutine real_ssygvd(a,b,w,uplo,jobz,itype)
581
582    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
583    !> overwritten on return anyway)
584    real(rsp), intent(inout) :: a(:,:)
585
586    !> contains the second matrix for the solver (overwritten by Cholesky factorization)
587    real(rsp), intent(inout) :: b(:,:)
588
589    !> eigenvalues
590    real(rsp), intent(out) :: w(:)
591
592    !> upper or lower triangle of the matrix
593    character, intent(in) :: uplo
594
595    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
596    character, intent(in) :: jobz
597
598    !> optional specifies the problem type to be solved 1:A*x=(lambda)*B*x, 2:A*B*x=(lambda)*x,
599    !> 3:B*A*x=(lambda)*x default is 1
600    integer, optional, intent(in) :: itype
601
602    real(rsp), allocatable :: work(:)
603    integer n, info, iitype
604    integer :: int_idealwork, iidealwork(1)
605    integer, allocatable :: iwork(:)
606    real(rsp) :: idealwork(1)
607    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
608    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
609    @:ASSERT(all(shape(a)==shape(b)))
610    @:ASSERT(all(shape(a)==size(w,dim=1)))
611    n=size(a,dim=1)
612    @:ASSERT(n>0)
613    if (present(itype)) then
614      iitype = itype
615    else
616      iitype = 1
617    end if
618    @:ASSERT(iitype >= 1 .and. iitype <= 3 )
619    call ssygvd(iitype, jobz, uplo, n, a, n, b, n, w, idealwork, -1, &
620         & iidealwork, -1, info)
621    if (info/=0) then
622       call error("Failure in SSYGVD to determine optimum workspace")
623    endif
624    int_idealwork=floor(idealwork(1))
625    allocate(work(int_idealwork))
626    allocate(iwork(iidealwork(1)))
627    call ssygvd(iitype, jobz, uplo, n, a, n, b, n, w, work, int_idealwork, &
628         & iwork, iidealwork(1), info)
629    if (info/=0) then
630       if (info<0) then
63199280 format ('Failure in diagonalisation routine ssygvd,', &
632          & ' illegal argument at position ',i6)
633          write(error_string, 99280) info
634          call error(error_string)
635       else if (info <= n) then
63699290 format ('Failure in diagonalisation routine ssygvd,', &
637          & ' diagonal element ',i6,' did not converge to zero.')
638          write(error_string, 99290) info
639          call error(error_string)
640       else
64199300 format ('Failure in diagonalisation routine ssygvd,', &
642          & ' non-positive definite overlap! Minor ',i6,' responsible.')
643          write(error_string, 99300) info - n
644          call error(error_string)
645       endif
646    endif
647
648  end subroutine real_ssygvd
649
650
651  !> Double precision eigensolver for generalized symmetric matrix problem divide and conquer
652  subroutine dble_dsygvd(a,b,w,uplo,jobz,itype)
653
654    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
655    !> overwritten on return anyway)
656    real(rdp), intent(inout) :: a(:,:)
657
658    !> contains the second matrix for the solver (overwritten by Cholesky factorization)
659    real(rdp), intent(inout) :: b(:,:)
660
661    !> eigenvalues
662    real(rdp), intent(out) :: w(:)
663
664    !> upper or lower triangle of the matrix
665    character, intent(in) :: uplo
666
667    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
668    character, intent(in) :: jobz
669
670    !> optional specifies the problem type to be solved 1:A*x=(lambda)*B*x, 2:A*B*x=(lambda)*x,
671    !> 3:B*A*x=(lambda)*x default is 1
672    integer, optional, intent(in) :: itype
673
674    real(rdp), allocatable :: work(:)
675    integer n, info, iitype
676    integer :: int_idealwork, iidealwork(1)
677    integer, allocatable :: iwork(:)
678    real(rdp) :: idealwork(1)
679    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
680    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
681    @:ASSERT(all(shape(a)==shape(b)))
682    @:ASSERT(all(shape(a)==size(w,dim=1)))
683    n=size(a,dim=1)
684    @:ASSERT(n>0)
685    if (present(itype)) then
686      iitype = itype
687    else
688      iitype = 1
689    end if
690    @:ASSERT(iitype >= 1 .and. iitype <= 3 )
691    call dsygvd(iitype, jobz, uplo, n, a, n, b, n, w, idealwork, -1, &
692        & iidealwork, -1, info)
693    if (info/=0) then
694      call error("Failure in DSYGVD to determine optimum workspace")
695    endif
696    int_idealwork=floor(idealwork(1))
697    allocate(work(int_idealwork))
698    allocate(iwork(iidealwork(1)))
699    call dsygvd(iitype, jobz, uplo, n, a, n, b, n, w, work, int_idealwork, &
700        & iwork, iidealwork(1), info)
701    if (info/=0) then
702      if (info<0) then
70399310 format ('Failure in diagonalisation routine dsygvd,', &
704          & ' illegal argument at position ',i6)
705        write(error_string, 99310) info
706        call error(error_string)
707      else if (info <= n) then
70899320 format ('Failure in diagonalisation routine dsygvd,', &
709          & ' diagonal element ',i6,' did not converge to zero.')
710        write(error_string, 99320) info
711        call error(error_string)
712      else
71399330 format ('Failure in diagonalisation routine dsygvd,', &
714          & ' non-positive definite overlap! Minor ',i6,' responsible.')
715        write(error_string, 99330) info - n
716        call error(error_string)
717      endif
718    endif
719
720  end subroutine dble_dsygvd
721
722
723  !> Complex eigensolver for generalized Hermitian matrix problem divide and conquer
724  subroutine cmplx_chegvd(a,b,w,uplo,jobz,itype)
725
726    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
727    !> overwritten on return anyway)
728    complex(rsp), intent(inout) :: a(:,:)
729
730    !> contains the second matrix for the solver (overwritten by Cholesky factorization)
731    complex(rsp), intent(inout) :: b(:,:)
732
733    !> eigenvalues
734    real(rsp), intent(out) :: w(:)
735
736    !> upper or lower triangle of the matrix
737    character, intent(in) :: uplo
738
739    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
740    character, intent(in) :: jobz
741
742    !> optional specifies the problem type to be solved 1:A*x=(lambda)*B*x, 2:A*B*x=(lambda)*x,
743    !> 3:B*A*x=(lambda)*x default is 1
744    integer, optional, intent(in) :: itype
745
746    complex(rsp), allocatable :: work(:)
747    real(rsp), allocatable :: rwork(:)
748    integer n, info, iitype
749    integer :: int_idealwork, int_ridealwork, iidealwork(1)
750    integer, allocatable :: iwork(:)
751    complex(rsp) :: idealwork(1)
752    real(rsp) :: ridealwork(1)
753    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
754    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
755    @:ASSERT(all(shape(a)==shape(b)))
756    @:ASSERT(all(shape(a)==size(w,dim=1)))
757    n=size(a,dim=1)
758    @:ASSERT(n>0)
759    if (present(itype)) then
760      iitype = itype
761    else
762      iitype = 1
763    end if
764    @:ASSERT(iitype >= 1 .and. iitype <= 3 )
765    call chegvd(iitype, jobz, uplo, n, a, n, b, n, w, idealwork, -1, &
766        & ridealwork, -1, iidealwork, -1, info)
767    if (info/=0) then
768      call error("Failure in CHEGVD to determine optimum workspace")
769    endif
770    int_idealwork=floor(real(idealwork(1)))
771    int_ridealwork=floor(ridealwork(1))
772    allocate(work(int_idealwork))
773    allocate(rwork(int_ridealwork))
774    allocate(iwork(iidealwork(1)))
775    call chegvd(iitype, jobz, uplo, n, a, n, b, n, w, work, int_idealwork,&
776        & rwork, int_ridealwork, iwork, iidealwork(1), info)
777    if (info/=0) then
778      if (info<0) then
77999340 format ('Failure in diagonalisation routine chegvd,', &
780          & ' illegal argument at position ',i6)
781        write(error_string, 99340) info
782        call error(error_string)
783      else if (info <= n) then
78499350 format ('Failure in diagonalisation routine chegvd,', &
785          & ' diagonal element ',i6,' did not converge to zero.')
786        write(error_string, 99350) info
787        call error(error_string)
788      else
78999360 format ('Failure in diagonalisation routine chegvd,', &
790          & ' non-positive definite overlap! Minor ',i6,' responsible.')
791        write(error_string, 99360) info - n
792        call error(error_string)
793      endif
794    endif
795
796  end subroutine cmplx_chegvd
797
798
799  !> Double complex eigensolver for generalized Hermitian matrix problem divide and conquer
800  subroutine dblecmplx_zhegvd(a,b,w,uplo,jobz,itype)
801
802    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
803    !> overwritten on return anyway)
804    complex(rdp), intent(inout) :: a(:,:)
805
806    !> contains the second matrix for the solver (overwritten by Cholesky factorization)
807    complex(rdp), intent(inout) :: b(:,:)
808
809    !> eigenvalues
810    real(rdp), intent(out) :: w(:)
811
812    !> upper or lower triangle of the matrix
813    character, intent(in) :: uplo
814
815    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
816    character, intent(in) :: jobz
817
818    !> optional specifies the problem type to be solved 1:A*x=(lambda)*B*x, 2:A*B*x=(lambda)*x,
819    !> 3:B*A*x=(lambda)*x default is 1
820    integer, optional, intent(in) :: itype
821
822    complex(rdp), allocatable :: work(:)
823    real(rdp), allocatable :: rwork(:)
824    integer n, info, iitype
825    integer :: int_idealwork, int_ridealwork, iidealwork(1)
826    integer, allocatable :: iwork(:)
827    complex(rdp) :: idealwork(1)
828    real(rdp) :: ridealwork(1)
829    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
830    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
831    @:ASSERT(all(shape(a)==shape(b)))
832    @:ASSERT(all(shape(a)==size(w,dim=1)))
833    n=size(a,dim=1)
834    @:ASSERT(n>0)
835    if (present(itype)) then
836      iitype = itype
837    else
838      iitype = 1
839    end if
840    @:ASSERT(iitype >= 1 .and. iitype <= 3 )
841    call zhegvd(iitype, jobz, uplo, n, a, n, b, n, w, idealwork, -1, &
842        & ridealwork, -1, iidealwork, -1, info)
843    if (info/=0) then
844      call error("Failure in ZHEGVD to determine optimum workspace")
845    endif
846    int_idealwork=floor(real(idealwork(1)))
847    int_ridealwork=floor(ridealwork(1))
848    allocate(work(int_idealwork))
849    allocate(rwork(int_ridealwork))
850    allocate(iwork(iidealwork(1)))
851    call zhegvd(iitype, jobz, uplo, n, a, n, b, n, w, work, int_idealwork, &
852        & rwork, int_ridealwork, iwork, iidealwork(1), info)
853    if (info/=0) then
854      if (info<0) then
85599370 format ('Failure in diagonalisation routine zhegvd,', &
856          & ' illegal argument at position ',i6)
857        write(error_string, 99370) info
858        call error(error_string)
859      else if (info <= n) then
86099380 format ('Failure in diagonalisation routine zhegvd,', &
861          & ' diagonal element ',i6,' did not converge to zero.')
862        write(error_string, 99380) info
863        call error(error_string)
864      else
86599390 format ('Failure in diagonalisation routine zhegvd,', &
866          & ' non-positive definite overlap! Minor ',i6,' responsible.')
867        write(error_string, 99390) info - n
868        call error(error_string)
869      endif
870    endif
871
872  end subroutine dblecmplx_zhegvd
873
874
875  !> Real eigensolver for generalized symmetric matrix problem - Relatively Robust
876  !> Representation, optionally use the subspace form if w is smaller than the size of a and b, then
877  !> only the first n eigenvalues/eigenvectors are found.
878  !> This version re-uses a triangle of a matrix (saving an additional allocation that was in the
879  !> previous version).
880  !> Based in part on deMon routine from T. Heine
881  subroutine real_ssygvr(a,b,w,uplo,jobz,itype,ilIn,iuIn)
882
883    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
884    !> overwritten on return anyway)
885    real(rsp), intent(inout) :: a(:,:)
886
887    !> contains the second matrix for the solver (overwritten by Cholesky factorization)
888    real(rsp), intent(inout) :: b(:,:)
889
890    !> eigenvalues
891    real(rsp), intent(out) :: w(:)
892
893    !> upper or lower triangle of both matrices
894    character, intent(in) :: uplo
895
896    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
897    character, intent(in) :: jobz
898
899    !> specifies the problem type to be solved 1:A*x=(lambda)*B*x, 2:A*B*x=(lambda)*x,
900    !> 3:B*A*x=(lambda)*x default is 1
901    integer, optional, intent(in) :: itype
902
903    !> lower range of eigenstates
904    integer, optional, intent(in) :: ilIn
905
906    !> upper range of eigenstates
907    integer, optional, intent(in) :: iuIn
908
909    real(rsp), allocatable :: work(:)
910    real(rsp) :: tmpWork(1)
911    real(rsp), allocatable :: tmpChole(:)
912    integer :: lwork
913    integer, allocatable :: iwork(:)
914    integer :: tmpIWork(1)
915    integer :: liwork
916    integer, allocatable :: isuppz(:)
917    integer :: n, info, iitype
918    integer :: m, neig
919    real(rsp) :: abstol
920    logical :: wantz,upper
921    character :: trans
922    real(rsp) :: vl, vu
923    integer :: il, iu
924    logical :: subspace
925    integer :: ii, jj
926    character :: uplo_new
927
928    n = size(a, dim=1)
929
930    @:ASSERT(n > 0)
931    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
932    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
933    @:ASSERT(all(shape(a)==shape(b)))
934    @:ASSERT(present(ilIn) .eqv. present(iuIn))
935
936    subspace = (size(w) < n)
937    if (subspace) then
938      if (present(ilIn)) then
939        @:ASSERT(ilIn <= iuIn)
940        @:ASSERT(ilIn > 0)
941        @:ASSERT(n >= iuIn)
942        @:ASSERT(size(w,dim=1) == (iuIn - ilIn + 1))
943        il = ilIn
944        iu = iuIn
945      else
946        il = 1
947        iu = size(w,dim=1)
948      end if
949    else
950      if (present(ilIn)) then
951        @:ASSERT(ilIn == 1 .and. iuIn == n)
952      end if
953      il = 1
954      iu = n
955    end if
956
957    if (present(itype)) then
958      iitype = itype
959    else
960      iitype = 1
961    end if
962    @:ASSERT(iitype >= 1 .and. iitype <= 3 )
963
964    allocate(isuppz(2*n))
965
966    allocate(tmpChole(size(a,dim=1)))
967
968    wantz = ( jobz == 'V' .or. jobz == 'v' )
969    upper = ( uplo == 'U' .or. uplo == 'u' )
970    abstol = slamch( 'Safe minimum' )
971
972    info = 0
973
974    if (subspace) then
975      call ssyevr(jobz,'I',uplo,n,a,size(a,dim=1),vl,vu,il,iu,abstol,m,&
976          & w,b,size(b,dim=1),isuppz,tmpWork,-1,tmpIWork,-1,info)
977    else
978      call ssyevr(jobz,'A',uplo,n,a,size(a,dim=1),vl,vu,il,iu,abstol,m,&
979          & w,b,size(b,dim=1),isuppz,tmpWork,-1,tmpIWork,-1,info)
980    end if
981
982    if (info/=0) then
983      call error("Failure in SSYGVR to determine optimum workspace")
984    endif
985    lwork = floor(tmpWork(1))
986    liwork = floor(real(tmpIWork(1)))
987    allocate(work(lwork))
988    allocate(iwork(liwork))
989
990    ! Form a Cholesky factorization of B.
991    call spotrf( uplo, n, b, n, info )
992    if( info /= 0 ) then
993      info = n + info
99499400 format ('Failure in diagonalisation routine SSYGVR,', &
995          & ' unable to complete Cholesky factorization of B ',i6)
996      write(error_string, 99400) info
997      call error(error_string)
998    end if
999    ! Transform problem to standard eigenvalue problem and solve.
1000    call ssygst( iitype, uplo, n, a, n, b, n, info )
1001
1002    if( info /= 0 ) then
1003      write(error_string, *)'Failure in SSYGVR to transform to standard',info
1004      call error(error_string)
1005    end if
1006
1007    if ( wantz ) then
1008      ! Save Cholesky factor in the other triangle of H and tmpChole
1009      do ii = 1, n
1010        tmpChole(ii) = b(ii,ii)
1011      end do
1012      if ( upper ) then
1013        do jj = 1, n
1014          do ii = jj+1, n
1015            a(ii,jj) = b(jj,ii)
1016          end do
1017        end do
1018      else
1019        do jj = 1, n
1020          do ii = 1, jj-1
1021            a(ii,jj) = b(jj,ii)
1022          end do
1023        end do
1024      end if
1025    end if
1026
1027    if (subspace) then
1028      call ssyevr( jobz, 'I', uplo, n, a, n, vl, vu, il, iu, &
1029          & abstol, m, w, b, n, isuppz, work, lwork, iwork, &
1030          & liwork, info )
1031    else
1032      call ssyevr( jobz, 'A', uplo, n, a, n, vl, vu, il, iu, &
1033          & abstol, m, w, b, n, isuppz, work, lwork, iwork, &
1034          & liwork, info )
1035    end if
1036    if( info /= 0 ) then
103799410 format ('Failure in SSYEVR ',I6)
1038      write(error_string, 99410) info
1039      call error(error_string)
1040    end if
1041
1042    if ( wantz ) then
1043      ! Backtransform eigenvectors to the original problem.
1044
1045      do ii = 1, n
1046        a(ii,ii) = tmpChole(ii)
1047      end do
1048
1049      if ( upper ) then
1050        uplo_new = 'L'
1051        upper = .false.
1052      else
1053        uplo_new = 'U'
1054        upper = .true.
1055      end if
1056
1057      neig = n
1058      if( info > 0 ) then
1059        neig = info - 1
1060      end if
1061      if( iitype == 1 .or. iitype == 2 ) then
1062        ! For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
1063        ! backtransform eigenvectors: x = inv(L)'*y or inv(U)*y          !'
1064        if( upper ) then
1065          trans = 'N'
1066        else
1067          trans = 'T'
1068        end if
1069        call strsm('Left',uplo_new,trans,'Non-unit',n,neig,1.0,A,n, B,n)
1070      else if( iitype == 3 ) then
1071        ! For B*A*x=(lambda)*x;
1072        ! backtransform eigenvectors: x = L*y or U'*y     !'
1073        if( upper ) then
1074          trans = 'T'
1075        else
1076          trans = 'N'
1077        end if
1078        call strmm('Left',uplo_new,trans,'Non-unit',n,neig,1.0,a,n, b,n)
1079      end if
1080      do ii = 1,m
1081        a( 1:n, ii) = b( 1:n, ii )
1082      end do
1083      a( 1:n, m+1:n ) = 0.0
1084    end if
1085
1086  end subroutine real_ssygvr
1087
1088
1089  !> Double precision eigensolver for generalized symmetric matrix problem - Relatively Robust
1090  !> Representation, optionally use the subspace form if w is smaller than the size of a and b, then
1091  !> only the first n eigenvalues/eigenvectors are found.
1092  !> This version re-uses a triangle of a matrix (saving an additional allocation that was in the
1093  !> previous version).
1094  !> Based in part on deMon routine from T. Heine
1095  subroutine dble_dsygvr(a,b,w,uplo,jobz,itype,ilIn,iuIn)
1096
1097    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
1098    !> overwritten on return anyway)
1099    real(rdp), intent(inout) :: a(:,:)
1100
1101    !> contains the second matrix for the solver (overwritten by Cholesky factorization)
1102    real(rdp), intent(inout) :: b(:,:)
1103
1104    !> eigenvalues
1105    real(rdp), intent(out) :: w(:)
1106
1107    !> upper or lower triangle of both matrices
1108    character, intent(in) :: uplo
1109
1110    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
1111    character, intent(in) :: jobz
1112
1113    !> specifies the problem type to be solved 1:A*x=(lambda)*B*x, 2:A*B*x=(lambda)*x,
1114    !> 3:B*A*x=(lambda)*x default is 1
1115    integer, optional, intent(in) :: itype
1116
1117    !> lower range of eigenstates
1118    integer, optional, intent(in) :: ilIn
1119
1120    !> upper range of eigenstates
1121    integer, optional, intent(in) :: iuIn
1122
1123    real(rdp), allocatable :: work(:)
1124    real(rdp) :: tmpWork(1)
1125    real(rdp), allocatable :: tmpChole(:)
1126    integer :: lwork
1127    integer, allocatable :: iwork(:)
1128    integer :: tmpIWork(1)
1129    integer :: liwork
1130    integer, allocatable :: isuppz(:)
1131    integer :: n, info, iitype
1132    integer :: m, neig
1133    real(rdp) :: abstol
1134    logical :: wantz,upper
1135    character :: trans
1136    real(rdp) :: vl, vu
1137    integer :: il, iu
1138    logical :: subspace
1139    integer :: ii, jj
1140    character :: uplo_new
1141
1142    n = size(a, dim=1)
1143
1144    @:ASSERT(n > 0)
1145    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
1146    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
1147    @:ASSERT(all(shape(a)==shape(b)))
1148    @:ASSERT(present(ilIn) .eqv. present(iuIn))
1149
1150    subspace = (size(w) < n)
1151    if (subspace) then
1152      if (present(ilIn)) then
1153        @:ASSERT(ilIn <= iuIn)
1154        @:ASSERT(ilIn > 0)
1155        @:ASSERT(n >= iuIn)
1156        @:ASSERT(size(w,dim=1) == (iuIn - ilIn + 1))
1157        il = ilIn
1158        iu = iuIn
1159      else
1160        il = 1
1161        iu = size(w,dim=1)
1162      end if
1163    else
1164      if (present(ilIn)) then
1165        @:ASSERT(ilIn == 1 .and. iuIn == n)
1166      end if
1167      il = 1
1168      iu = n
1169    end if
1170
1171    if (present(itype)) then
1172      iitype = itype
1173    else
1174      iitype = 1
1175    end if
1176    @:ASSERT(iitype >= 1 .and. iitype <= 3 )
1177
1178    allocate(isuppz(2*n))
1179
1180    allocate(tmpChole(size(a,dim=1)))
1181
1182    wantz = ( jobz == 'V' .or. jobz == 'v' )
1183    upper = ( uplo == 'U' .or. uplo == 'u' )
1184    abstol = dlamch( 'Safe minimum' )
1185
1186    info = 0
1187
1188    if (subspace) then
1189      call dsyevr(jobz,'I',uplo,n,a,size(a,dim=1),vl,vu,il,iu,abstol,m,&
1190          & w,b,size(b,dim=1),isuppz,tmpWork,-1,tmpIWork,-1,info)
1191    else
1192      call dsyevr(jobz,'A',uplo,n,a,size(a,dim=1),vl,vu,il,iu,abstol,m,&
1193          & w,b,size(b,dim=1),isuppz,tmpWork,-1,tmpIWork,-1,info)
1194    end if
1195
1196    if (info/=0) then
1197      call error("Failure in DSYGVR to determine optimum workspace")
1198    endif
1199    lwork = floor(tmpWork(1))
1200    liwork = floor(real(tmpIWork(1)))
1201    allocate(work(lwork))
1202    allocate(iwork(liwork))
1203
1204    ! Form a Cholesky factorization of B.
1205    call dpotrf( uplo, n, b, n, info )
1206    if( info /= 0 ) then
1207      info = n + info
120899400 format ('Failure in diagonalisation routine DSYGVR,', &
1209          & ' unable to complete Cholesky factorization of B ',i6)
1210      write(error_string, 99400) info
1211      call error(error_string)
1212    end if
1213    ! Transform problem to standard eigenvalue problem and solve.
1214    call dsygst( iitype, uplo, n, a, n, b, n, info )
1215
1216    if( info /= 0 ) then
1217      write(error_string, *)'Failure in DSYGVR to transform to standard',info
1218      call error(error_string)
1219    end if
1220
1221    if ( wantz ) then
1222      ! Save Cholesky factor in the other triangle of H and tmpChole
1223      do ii = 1, n
1224        tmpChole(ii) = b(ii,ii)
1225      end do
1226      if ( upper ) then
1227        do jj = 1, n
1228          do ii = jj + 1, n
1229            a(ii,jj) = b(jj,ii)
1230          end do
1231        end do
1232      else
1233        do jj = 1, n
1234          do ii = 1, jj - 1
1235            a( ii, jj ) = b( jj, ii )
1236          end do
1237        end do
1238      end if
1239    end if
1240
1241    if (subspace) then
1242      call dsyevr( jobz, 'I', uplo, n, a, n, vl, vu, il, iu, &
1243          & abstol, m, w, b, n, isuppz, work, lwork, iwork, &
1244          & liwork, info )
1245    else
1246      call dsyevr( jobz, 'A', uplo, n, a, n, vl, vu, il, iu, &
1247          & abstol, m, w, b, n, isuppz, work, lwork, iwork, &
1248          & liwork, info )
1249    end if
1250    if( info /= 0 ) then
125199410 format ('Failure in DSYEVR ',I6)
1252      write(error_string, 99410) info
1253      call error(error_string)
1254    end if
1255
1256    if ( wantz ) then
1257      ! Backtransform eigenvectors to the original problem.
1258
1259      do ii = 1, n
1260        a(ii,ii) = tmpChole(ii)
1261      end do
1262
1263      if ( upper ) then
1264        uplo_new = 'L'
1265        upper = .false.
1266      else
1267        uplo_new = 'U'
1268        upper = .true.
1269      end if
1270
1271      neig = n
1272      if( info > 0 ) then
1273        neig = info - 1
1274      end if
1275      if( iitype == 1 .or. iitype == 2 ) then
1276        ! For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
1277        ! backtransform eigenvectors: x = inv(L)'*y or inv(U)*y          !'
1278        if( upper ) then
1279          trans = 'N'
1280        else
1281          trans = 'T'
1282        end if
1283        call dtrsm('Left',uplo_new,trans,'Non-unit',n,neig,1.0d0,A,n, B,n)
1284      else if( iitype == 3 ) then
1285        ! For B*A*x=(lambda)*x;
1286        ! backtransform eigenvectors: x = L*y or U'*y               !'
1287        if( upper ) then
1288          trans = 'T'
1289        else
1290          trans = 'N'
1291        end if
1292        call dtrmm('Left',uplo_new,trans,'Non-unit',n,neig,1.0d0,a,n, b,n)
1293      end if
1294      do ii = 1,m
1295        a( 1:n, ii) = b( 1:n, ii )
1296      end do
1297      a( 1:n, m+1:n ) = 0.0d0
1298    end if
1299
1300  end subroutine dble_dsygvr
1301
1302
1303  !> Complex precision eigensolver for generalized symmetric matrix problem - Relatively Robust
1304  !> Representation, optionally use the subspace form if w is smaller than the size of a and b, then
1305  !> only the first n eigenvalues/eigenvectors are found.
1306  !> This version re-uses a triangle of a matrix (saving an additional allocation that was in the
1307  !> previous version).
1308  !> Based in part on deMon routine from T. Heine
1309  subroutine cmplx_chegvr(a,b,w,uplo,jobz,itype,ilIn,iuIn)
1310
1311    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
1312    !> overwritten on return anyway)
1313    complex(rsp), intent(inout) :: a(:,:)
1314
1315    !> contains the second matrix for the solver (overwritten by Cholesky factorization)
1316    complex(rsp), intent(inout) :: b(:,:)
1317
1318    !> eigenvalues
1319    real(rsp), intent(out) :: w(:)
1320
1321    !> upper or lower triangle of both matrices
1322    character, intent(in) :: uplo
1323
1324    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
1325    character, intent(in) :: jobz
1326
1327    !> specifies the problem type to be solved 1:A*x=(lambda)*B*x, 2:A*B*x=(lambda)*x,
1328    !> 3:B*A*x=(lambda)*x default is 1
1329    integer, optional, intent(in) :: itype
1330
1331    !> lower range of eigenstates
1332    integer, optional, intent(in) :: ilIn
1333
1334    !> upper range of eigenstates
1335    integer, optional, intent(in) :: iuIn
1336
1337    complex(rsp), allocatable :: work(:)
1338    complex(rsp) :: tmpWork(1)
1339    real(rsp), allocatable :: rWork(:)
1340    real(rsp) :: tmpRWork(1)
1341    complex(rsp), allocatable :: tmpChole(:)
1342    integer :: lwork
1343    integer, allocatable :: iwork(:)
1344    integer :: tmpIWork(1)
1345    integer :: liwork
1346    integer :: lrwork
1347    integer, allocatable :: isuppz(:)
1348    integer :: n, info, iitype
1349    integer :: m, neig
1350    real(rsp) :: abstol
1351    logical :: wantz,upper
1352    character :: trans
1353    real(rsp) :: vl, vu
1354    integer :: il, iu
1355    logical :: subspace
1356    integer :: ii, jj
1357    character :: uplo_new
1358
1359    n = size(a, dim=1)
1360    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
1361    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
1362    @:ASSERT(all(shape(a)==shape(b)))
1363
1364    @:ASSERT(present(ilIn) .eqv. present(iuIn))
1365
1366    subspace = (size(w) < n)
1367    if (subspace) then
1368      if (present(ilIn)) then
1369        @:ASSERT(ilIn <= iuIn)
1370        @:ASSERT(ilIn > 0)
1371        @:ASSERT(n >= iuIn)
1372        @:ASSERT(size(w,dim=1) == (iuIn - ilIn + 1))
1373        il = ilIn
1374        iu = iuIn
1375      else
1376        il = 1
1377        iu = size(w,dim=1)
1378      end if
1379    else
1380      if (present(ilIn)) then
1381        @:ASSERT(ilIn == 1 .and. iuIn == n)
1382      end if
1383      il = 1
1384      iu = n
1385    end if
1386
1387    if (present(itype)) then
1388      iitype = itype
1389    else
1390      iitype = 1
1391    end if
1392    @:ASSERT(iitype >= 1 .and. iitype <= 3 )
1393
1394    allocate(isuppz(2*n))
1395
1396    allocate(tmpChole(size(a,dim=1)))
1397
1398    wantz = ( jobz == 'V' .or. jobz == 'v' )
1399    upper = ( uplo == 'U' .or. uplo == 'u' )
1400    abstol = SLAMCH( 'Safe minimum' )
1401
1402    info = 0
1403
1404    if (subspace) then
1405      call cheevr(jobz,'I',uplo,n,a,size(a,dim=1),vl,vu,il,iu,abstol,m,&
1406          & w,b,size(b,dim=1),isuppz,tmpWork,-1,tmpRwork,-1,tmpIWork,-1,info)
1407    else
1408      call cheevr(jobz,'A',uplo,n,a,size(a,dim=1),vl,vu,il,iu,abstol,m,&
1409          & w,b,size(b,dim=1),isuppz,tmpWork,-1,tmpRwork,-1,tmpIWork,-1,info)
1410    end if
1411
1412    if (info/=0) then
1413      call error("Failure in CHEEVR to determine optimum workspace")
1414    endif
1415    lwork = floor(real(tmpWork(1)))
1416    liwork = floor(real(tmpIWork(1)))
1417    lrwork = floor(tmpRwork(1))
1418    allocate(work(lwork))
1419    allocate(iwork(liwork))
1420    allocate(rwork(lrwork))
1421
1422    ! Form a Cholesky factorization of B.
1423    call cpotrf( uplo, n, b, n, info )
1424    if( info /= 0 ) then
1425      info = n + info
142699400 format ('Failure in diagonalisation routine CHEGVR,', &
1427          & ' unable to complete Cholesky factorization of B ',i6)
1428      write(error_string, 99400) info
1429      call error(error_string)
1430    end if
1431    ! Transform problem to standard eigenvalue problem and solve.
1432    call chegst( iitype, uplo, n, a, n, b, n, info )
1433
1434    if( info /= 0 ) then
1435      write(error_string, *)'Failure in CHEGST to transform to standard',info
1436      call error(error_string)
1437    end if
1438
1439    if ( wantz ) then
1440      ! Save Cholesky factor in the other triangle of H and tmpChole
1441      do ii = 1, n
1442        tmpChole(ii) = b(ii,ii)
1443      end do
1444      if ( upper ) then
1445        do jj = 1, n
1446          do ii = jj+1, n
1447            a(ii,jj) = conjg(b(jj,ii))
1448          end do
1449        end do
1450      else
1451        do jj = 1, n
1452          do ii = 1, jj-1
1453            a(ii,jj) = conjg(b(jj,ii))
1454          end do
1455        end do
1456      end if
1457    end if
1458
1459    if (subspace) then
1460      call cheevr( jobz, 'I', uplo, n, a, n, vl, vu, il, iu, &
1461          & abstol, m, w, b, n, isuppz, work, lwork, rwork, lrwork, iwork, &
1462          & liwork, info )
1463    else
1464      call cheevr( jobz, 'A', uplo, n, a, n, vl, vu, il, iu, &
1465          & abstol, m, w, b, n, isuppz, work, lwork, rwork, lrwork, iwork, &
1466          & liwork, info )
1467    end if
1468    if( info /= 0 ) then
146999410 format ('Failure in CHEEVR ',I6)
1470      write(error_string, 99410) info
1471      call error(error_string)
1472    end if
1473
1474    if ( wantz ) then
1475      ! Backtransform eigenvectors to the original problem.
1476
1477      do ii = 1, n
1478        a(ii,ii) = tmpChole(ii)
1479      end do
1480
1481      if ( upper ) then
1482        uplo_new = 'L'
1483        upper = .false.
1484      else
1485        uplo_new = 'U'
1486        upper = .true.
1487      end if
1488
1489      neig = n
1490      if( info > 0 ) then
1491        neig = info - 1
1492      end if
1493      if( iitype == 1 .or. iitype == 2 ) then
1494        ! For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
1495        ! backtransform eigenvectors: x = inv(L)'*y or inv(U)*y          !'
1496        if( upper ) then
1497          trans = 'N'
1498        else
1499          trans = 'C'
1500        end if
1501        call ctrsm('Left',uplo_new,trans,'Non-unit',n,neig, &
1502            & cmplx(1.0,0.0,rsp),A,n, B,n)
1503      else if( iitype == 3 ) then
1504        ! For B*A*x=(lambda)*x;
1505        ! backtransform eigenvectors: x = L*y or U'*y               !'
1506        if( upper ) then
1507          trans = 'C'
1508        else
1509          trans = 'N'
1510        end if
1511        call ctrmm('Left',uplo_new,trans,'Non-unit',n,neig, &
1512            & cmplx(1.0,0.0,rsp),a,n, b,n)
1513      end if
1514      do ii = 1,m
1515        a( 1:n, ii) = b( 1:n, ii )
1516      end do
1517      a( 1:n, m+1:n ) = cmplx(0.0,0.0,rsp)
1518    end if
1519
1520  end subroutine cmplx_chegvr
1521
1522
1523  !> Double complex precision eigensolver for generalized symmetric matrix problem - Relatively
1524  !> Robust Representation, optionally use the subspace form if w is smaller than the size of a and
1525  !> b, then only the first n eigenvalues/eigenvectors are found.
1526  !> This version re-uses a triangle of a matrix (saving an additional allocation that was in the
1527  !> previous version).
1528  !> Based in part on deMon routine from T. Heine
1529  subroutine dblecmplx_zhegvr(a,b,w,uplo,jobz,itype,ilIn,iuIn)
1530
1531    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
1532    !> overwritten on return anyway)
1533    complex(rdp), intent(inout) :: a(:,:)
1534
1535    !> contains the second matrix for the solver (overwritten by Cholesky factorization)
1536    complex(rdp), intent(inout) :: b(:,:)
1537
1538    !> eigenvalues
1539    real(rdp), intent(out) :: w(:)
1540
1541    !> upper or lower triangle of both matrices
1542    character, intent(in) :: uplo
1543
1544    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
1545    character, intent(in) :: jobz
1546
1547    !> specifies the problem type to be solved 1:A*x=(lambda)*B*x, 2:A*B*x=(lambda)*x,
1548    !> 3:B*A*x=(lambda)*x default is 1
1549    integer, optional, intent(in) :: itype
1550
1551    !> lower range of eigenstates
1552    integer, optional, intent(in) :: ilIn
1553
1554    !> upper range of eigenstates
1555    integer, optional, intent(in) :: iuIn
1556
1557    complex(rdp), allocatable :: work(:)
1558    complex(rdp) :: tmpWork(1)
1559    real(rdp), allocatable :: rWork(:)
1560    real(rdp) :: tmpRWork(1)
1561    complex(rdp), allocatable :: tmpChole(:)
1562    integer :: lwork
1563    integer, allocatable :: iwork(:)
1564    integer :: tmpIWork(1)
1565    integer :: liwork
1566    integer :: lrwork
1567    integer, allocatable :: isuppz(:)
1568    integer :: n, info, iitype
1569    integer :: m, neig
1570    real(rdp) :: abstol
1571    logical :: wantz,upper
1572    character :: trans
1573    real(rdp) :: vl, vu
1574    integer :: il, iu
1575    logical :: subspace
1576    integer :: ii, jj
1577    character :: uplo_new
1578
1579    n = size(a, dim=1)
1580
1581    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
1582    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
1583    @:ASSERT(all(shape(a)==shape(b)))
1584    @:ASSERT(present(ilIn) .eqv. present(iuIn))
1585
1586    subspace = (size(w) < n)
1587    if (subspace) then
1588      if (present(ilIn)) then
1589        @:ASSERT(ilIn <= iuIn)
1590        @:ASSERT(ilIn > 0)
1591        @:ASSERT(n >= iuIn)
1592        @:ASSERT(size(w,dim=1) == (iuIn - ilIn + 1))
1593        il = ilIn
1594        iu = iuIn
1595      else
1596        il = 1
1597        iu = size(w,dim=1)
1598      end if
1599    else
1600      if (present(ilIn)) then
1601        @:ASSERT(ilIn == 1 .and. iuIn == n)
1602      end if
1603      il = 1
1604      iu = n
1605    end if
1606
1607    if (present(itype)) then
1608      iitype = itype
1609    else
1610      iitype = 1
1611    end if
1612    @:ASSERT(iitype >= 1 .and. iitype <= 3 )
1613
1614    allocate(isuppz(2*n))
1615
1616    allocate(tmpChole(size(a,dim=1)))
1617
1618    wantz = ( jobz == 'V' .or. jobz == 'v' )
1619    upper = ( uplo == 'U' .or. uplo == 'u' )
1620    abstol = DLAMCH( 'Safe minimum' )
1621
1622    info = 0
1623
1624    if (subspace) then
1625      call zheevr(jobz,'I',uplo,n,a,size(a,dim=1),vl,vu,il,iu,abstol,m,&
1626          & w,b,size(b,dim=1),isuppz,tmpWork,-1,tmpRwork,-1,tmpIWork,-1,info)
1627    else
1628      call zheevr(jobz,'A',uplo,n,a,size(a,dim=1),vl,vu,il,iu,abstol,m,&
1629          & w,b,size(b,dim=1),isuppz,tmpWork,-1,tmpRwork,-1,tmpIWork,-1,info)
1630    end if
1631
1632    if (info/=0) then
1633      call error("Failure in ZHEEVR to determine optimum workspace")
1634    endif
1635    lwork = floor(real(tmpWork(1)))
1636    liwork = floor(real(tmpIWork(1)))
1637    lrwork = floor(tmpRwork(1))
1638    allocate(work(lwork))
1639    allocate(iwork(liwork))
1640    allocate(rwork(lrwork))
1641
1642    ! Form a Cholesky factorization of B.
1643    call zpotrf( uplo, n, b, n, info )
1644    if( info /= 0 ) then
1645      info = n + info
164699400 format ('Failure in diagonalisation routine ZHEEVR,', &
1647          & ' unable to complete Cholesky factorization of B ',i6)
1648      write(error_string, 99400) info
1649      call error(error_string)
1650    end if
1651    ! Transform problem to standard eigenvalue problem and solve.
1652    call zhegst( iitype, uplo, n, a, n, b, n, info )
1653
1654    if( info /= 0 ) then
1655      write(error_string, *)'Failure in ZHEGST to transform to standard',info
1656      call error(error_string)
1657    end if
1658
1659    if ( wantz ) then
1660      ! Save Cholesky factor in the other triangle of H and tmpChole
1661      do ii = 1, n
1662        tmpChole(ii) = b(ii,ii)
1663      end do
1664      if ( upper ) then
1665        do jj = 1, n
1666          do ii = jj+1, n
1667            a(ii,jj) = conjg(b(jj,ii))
1668          end do
1669        end do
1670      else
1671        do jj = 1, n
1672          do ii = 1, jj-1
1673            a(ii,jj) = conjg(b(jj,ii))
1674          end do
1675        end do
1676      end if
1677    end if
1678
1679    if (subspace) then
1680      call zheevr( jobz, 'I', uplo, n, a, n, vl, vu, il, iu, &
1681          & abstol, m, w, b, n, isuppz, work, lwork, rwork, lrwork, iwork, &
1682          & liwork, info )
1683    else
1684      call zheevr( jobz, 'A', uplo, n, a, n, vl, vu, il, iu, &
1685          & abstol, m, w, b, n, isuppz, work, lwork, rwork, lrwork, iwork, &
1686          & liwork, info )
1687    end if
1688    if( info /= 0 ) then
168999410 format ('Failure in ZHEEVR ',I6)
1690      write(error_string, 99410) info
1691      call error(error_string)
1692    end if
1693
1694    if ( wantz ) then
1695      ! Backtransform eigenvectors to the original problem.
1696
1697      do ii = 1, n
1698        a(ii,ii) = tmpChole(ii)
1699      end do
1700
1701      if ( upper ) then
1702        uplo_new = 'L'
1703        upper = .false.
1704      else
1705        uplo_new = 'U'
1706        upper = .true.
1707      end if
1708
1709      neig = n
1710      if( info > 0 ) then
1711        neig = info - 1
1712      end if
1713      if( iitype == 1 .or. iitype == 2 ) then
1714        ! For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
1715        ! backtransform eigenvectors: x = inv(L)'*y or inv(U)*y          !'
1716        if( upper ) then
1717          trans = 'N'
1718        else
1719          trans = 'C'
1720        end if
1721        call ztrsm('Left',uplo_new,trans,'Non-unit',n,neig, &
1722            & cmplx(1.0,0.0,rdp),A,n, B,n)
1723      else if( iitype == 3 ) then
1724        ! For B*A*x=(lambda)*x;
1725        ! backtransform eigenvectors: x = L*y or U'*y               !'
1726        if( upper ) then
1727          trans = 'C'
1728        else
1729          trans = 'N'
1730        end if
1731        call ztrmm('Left',uplo_new,trans,'Non-unit',n,neig, &
1732            & cmplx(1.0,0.0,rdp),a,n, b,n)
1733      end if
1734      do ii = 1,m
1735        a( 1:n, ii) = b( 1:n, ii )
1736      end do
1737      a( 1:n, m+1:n ) = cmplx(0.0,0.0,rdp)
1738    end if
1739
1740  end subroutine dblecmplx_zhegvr
1741
1742
1743  !> Single precision banded symmetric generalised matrix eigensolver
1744  subroutine real_ssbgv(ab, bb, w, uplo, z)
1745
1746    !> contains the matrix for the solver (overwritten before exit)
1747    real(rsp), intent(inout) :: ab(:,:)
1748
1749    !> contains the second matrix for the solver (overwritten by split Cholesky factorization)
1750    real(rsp), intent(inout) :: bb(:,:)
1751
1752    !> eigenvalues
1753    real(rsp), intent(out) :: w(:)
1754
1755    !> upper or lower triangle of both matrices
1756    character, intent(in) :: uplo
1757
1758    !> returns calculated eigenvectors if present
1759    real(rsp), optional, intent(out) :: z(:,:)
1760
1761    real(rsp), allocatable :: work(:)
1762    integer :: n, ka, kb, ldab, ldbb, ldz, info
1763    character :: jobz
1764    real(rsp) :: zTmp(1,1)
1765
1766    if (present(z)) then
1767      jobz = 'v'
1768    else
1769      jobz = 'n'
1770    end if
1771
1772    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
1773    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
1774    @:ASSERT(all(shape(ab)==shape(bb)))
1775
1776    n = size(ab,dim=2)
1777
1778    ldab = size(ab,dim=1)
1779    ldbb = size(bb,dim=1)
1780    ka = ldab - 1
1781    kb = ldbb - 1
1782    @:ASSERT(ka >= 0)
1783    @:ASSERT(kb >= 0)
1784
1785    if (present(z)) then
1786      ldz = n
1787    else
1788      ldz = 1
1789    end if
1790
1791    allocate(work(3*n))
1792
1793    if (present(z)) then
1794      call ssbgv( jobz,uplo,n,ka,kb,ab,ldab,bb,ldbb,w,z,ldz,work,info )
1795    else
1796      call ssbgv( jobz,uplo,n,ka,kb,ab,ldab,bb,ldbb,w,zTmp,ldz,work,info )
1797    end if
1798
1799    if (info/=0) then
1800      if (info<0) then
180199440 format ('Failure in diagonalisation routine ssbgv,', &
1802          & ' illegal argument at position ',i6)
1803          write(error_string, 99440) info
1804          call error(error_string)
1805       else if (info <= n) then
180699450 format ('Failure in diagonalisation routine ssbgv,', &
1807          & ' tri diagonal element ',i6,' did not converge to zero.')
1808          write(error_string, 99450) info
1809          call error(error_string)
1810       else
181199460 format ('Failure in diagonalisation routine ssbgv,', &
1812          & ' non-positive definite overlap! ',i6)
1813          write(error_string, 99460) info
1814          call error(error_string)
1815       endif
1816    endif
1817
1818  end subroutine real_ssbgv
1819
1820
1821  !> Double precision banded symmetric generalised matrix eigensolver
1822  subroutine dble_dsbgv(ab, bb, w, uplo, z)
1823
1824    !> contains the matrix for the solver (overwritten before exit)
1825    real(rdp), intent(inout) :: ab(:,:)
1826
1827    !> contains the second matrix for the solver (overwritten by split Cholesky factorization)
1828    real(rdp), intent(inout) :: bb(:,:)
1829
1830    !> eigenvalues
1831    real(rdp), intent(out) :: w(:)
1832
1833    !> upper or lower triangle of both matrices
1834    character, intent(in) :: uplo
1835
1836    !> returns calculated eigenvectors if present
1837    real(rdp), optional, intent(out) :: z(:,:)
1838
1839    real(rdp), allocatable :: work(:)
1840    integer :: n, ka, kb, ldab, ldbb, ldz, info
1841    character :: jobz
1842    real(rdp) :: zTmp(1,1)
1843
1844    if (present(z)) then
1845      jobz = 'v'
1846    else
1847      jobz = 'n'
1848    end if
1849
1850    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
1851    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
1852    @:ASSERT(all(shape(ab)==shape(bb)))
1853
1854    n = size(ab,dim=2)
1855
1856    ldab = size(ab,dim=1)
1857    ldbb = size(bb,dim=1)
1858    ka = ldab - 1
1859    kb = ldbb - 1
1860    @:ASSERT(ka >= 0)
1861    @:ASSERT(kb >= 0)
1862
1863    if (present(z)) then
1864      ldz = n
1865    else
1866      ldz = 1
1867    end if
1868
1869    allocate(work(3*n))
1870
1871    if (present(z)) then
1872      call dsbgv( jobz,uplo,n,ka,kb,ab,ldab,bb,ldbb,w,z,ldz,work,info )
1873    else
1874      call dsbgv( jobz,uplo,n,ka,kb,ab,ldab,bb,ldbb,w,zTmp,ldz,work,info )
1875    end if
1876
1877    if (info/=0) then
1878      if (info<0) then
187999470 format ('Failure in diagonalisation routine dsbgv,', &
1880          & ' illegal argument at position ',i6)
1881          write(error_string, 99470) info
1882          call error(error_string)
1883       else if (info <= n) then
188499480 format ('Failure in diagonalisation routine dsbgv,', &
1885          & ' tri diagonal element ',i6,' did not converge to zero.')
1886          write(error_string, 99480) info
1887          call error(error_string)
1888       else
188999490 format ('Failure in diagonalisation routine dsbgv,', &
1890          & ' non-positive definite overlap! ',i6)
1891          write(error_string, 99490) info
1892          call error(error_string)
1893       endif
1894    endif
1895
1896  end subroutine dble_dsbgv
1897
1898
1899  !> Complex banded symmetric generalised matrix eigensolver
1900  subroutine cmplx_chbgv(ab, bb, w, uplo, z)
1901
1902    !> contains the matrix for the solver (overwritten before exit)
1903    complex(rsp), intent(inout) :: ab(:,:)
1904
1905    !> contains the second matrix for the solver (overwritten by split Cholesky factorization)
1906    complex(rsp), intent(inout) :: bb(:,:)
1907
1908    !> eigenvalues
1909    real(rsp), intent(out) :: w(:)
1910
1911    !> upper or lower triangle of both matrices
1912    character, intent(in) :: uplo
1913
1914    !> returns calculated eigenvectors if present
1915    complex(rsp), optional, intent(out) :: z(:,:)
1916
1917    complex(rsp), allocatable :: work(:)
1918    real(rsp), allocatable :: rwork(:)
1919    integer :: n, ka, kb, ldab, ldbb, ldz, info
1920    character :: jobz
1921    complex(rsp) :: zTmp(1,1)
1922
1923    if (present(z)) then
1924      jobz = 'v'
1925    else
1926      jobz = 'n'
1927    end if
1928
1929    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
1930    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
1931    @:ASSERT(all(shape(ab)==shape(bb)))
1932
1933    n = size(ab,dim=2)
1934
1935    ldab = size(ab,dim=1)
1936    ldbb = size(bb,dim=1)
1937    ka = ldab - 1
1938    kb = ldbb - 1
1939    @:ASSERT(ka >= 0)
1940    @:ASSERT(kb >= 0)
1941
1942    if (present(z)) then
1943      ldz = n
1944    else
1945      ldz = 1
1946    end if
1947
1948    allocate(work(n))
1949    allocate(rwork(3*n))
1950
1951    if (present(z)) then
1952      call chbgv( jobz,uplo,n,ka,kb,ab,ldab,bb,ldbb,w,z,ldz,work,rwork,info )
1953    else
1954      call chbgv( jobz,uplo,n,ka,kb,ab,ldab,bb,ldbb,w,zTmp,ldz,work,rwork, &
1955          & info )
1956    end if
1957
1958    if (info/=0) then
1959      if (info<0) then
196099500 format ('Failure in diagonalisation routine chbgv,', &
1961          & ' illegal argument at position ',i6)
1962          write(error_string, 99500) info
1963          call error(error_string)
1964       else if (info <= n) then
196599510 format ('Failure in diagonalisation routine chbgv,', &
1966          & ' tri diagonal element ',i6,' did not converge to zero.')
1967          write(error_string, 99510) info
1968          call error(error_string)
1969       else
197099520 format ('Failure in diagonalisation routine chbgv,', &
1971          & ' non-positive definite overlap! ',i6)
1972          write(error_string, 99520) info
1973          call error(error_string)
1974       endif
1975    endif
1976
1977  end subroutine cmplx_chbgv
1978
1979
1980  !> Complex double precision banded symmetric generalised matrix eigensolver
1981  subroutine dblecmplx_zhbgv(ab, bb, w, uplo, z)
1982
1983    !> contains the matrix for the solver (overwritten before exit)
1984    complex(rdp), intent(inout) :: ab(:,:)
1985
1986    !> contains the second matrix for the solver (overwritten by split Cholesky factorization)
1987    complex(rdp), intent(inout) :: bb(:,:)
1988
1989    !> eigenvalues
1990    real(rdp), intent(out) :: w(:)
1991
1992    !> upper or lower triangle of both matrices
1993    character, intent(in) :: uplo
1994
1995    !> returns calculated eigenvectors if present
1996    complex(rdp), optional, intent(out) :: z(:,:)
1997
1998    complex(rdp), allocatable :: work(:)
1999    real(rdp), allocatable :: rwork(:)
2000    integer :: n, ka, kb, ldab, ldbb, ldz, info
2001    character :: jobz
2002    complex(rdp) :: zTmp(1,1)
2003
2004    if (present(z)) then
2005      jobz = 'v'
2006    else
2007      jobz = 'n'
2008    end if
2009
2010    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
2011    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
2012    @:ASSERT(all(shape(ab)==shape(bb)))
2013
2014    n = size(ab,dim=2)
2015
2016    ldab = size(ab,dim=1)
2017    ldbb = size(bb,dim=1)
2018    ka = ldab - 1
2019    kb = ldbb - 1
2020    @:ASSERT(ka >= 0)
2021    @:ASSERT(kb >= 0)
2022
2023    if (present(z)) then
2024      ldz = n
2025    else
2026      ldz = 1
2027    end if
2028
2029    allocate(work(n))
2030    allocate(rwork(3*n))
2031
2032    if (present(z)) then
2033      call zhbgv( jobz,uplo,n,ka,kb,ab,ldab,bb,ldbb,w,z,ldz,work,rwork,info )
2034    else
2035      call zhbgv( jobz,uplo,n,ka,kb,ab,ldab,bb,ldbb,w,zTmp,ldz,work,rwork, &
2036          & info )
2037    end if
2038
2039    if (info/=0) then
2040      if (info<0) then
204199530 format ('Failure in diagonalisation routine zhbgv,', &
2042          & ' illegal argument at position ',i6)
2043          write(error_string, 99530) info
2044          call error(error_string)
2045       else if (info <= n) then
204699540 format ('Failure in diagonalisation routine zhbgv,', &
2047          & ' tri diagonal element ',i6,' did not converge to zero.')
2048          write(error_string, 99540) info
2049          call error(error_string)
2050       else
205199550 format ('Failure in diagonalisation routine zhbgv,', &
2052          & ' non-positive definite overlap! ',i6)
2053          write(error_string, 99550) info
2054          call error(error_string)
2055       endif
2056    endif
2057
2058  end subroutine dblecmplx_zhbgv
2059
2060
2061#:if WITH_GPU
2062
2063#: for DTYPE, VPREC, VTYPE, NAME in [('real', 's', 'real', 'ssygvd'),&
2064  & ('dble', 'd', 'real', 'dsygvd'), ('cmplx', 's', 'complex', 'chegvd'),&
2065  & ('dblecmplx', 'd', 'complex', 'zhegvd')]
2066  !> Generalised eigensolution for symmetric/hermitian matrices on GPU(s)
2067  subroutine ${DTYPE}$_magma_${NAME}$(ngpus, a, b, w, uplo, jobz, itype)
2068
2069    !> Number of GPUs to use
2070    integer, intent(in) :: ngpus
2071
2072    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
2073    !> overwritten on return anyway)
2074    ${VTYPE}$(r${VPREC}$p), intent(inout) :: a(:,:)
2075
2076    !> contains the second matrix for the solver (overwritten by Cholesky factorization)
2077    ${VTYPE}$(r${VPREC}$p), intent(inout) :: b(:,:)
2078
2079    !> eigenvalues
2080    real(r${VPREC}$p), intent(out) :: w(:)
2081
2082    !> upper or lower triangle of the matrix
2083    character, intent(in) :: uplo
2084
2085    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
2086    character, intent(in) :: jobz
2087
2088    !> optional specifies the problem type to be solved 1:A*x=(lambda)*B*x, 2:A*B*x=(lambda)*x,
2089    !> 3:B*A*x=(lambda)*x default is 1
2090    integer, optional, intent(in) :: itype
2091
2092    ${VTYPE}$(r${VPREC}$p), allocatable :: work(:)
2093
2094  #:if VTYPE == 'complex'
2095    real(r${VPREC}$p), allocatable :: rwork(:)
2096    integer :: lrwork
2097  #:endif
2098
2099    integer, allocatable :: iwork(:)
2100    integer :: lwork, liwork, n, info, iitype
2101
2102    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
2103    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
2104    @:ASSERT(all(shape(a)==shape(b)))
2105    @:ASSERT(all(shape(a)==size(w,dim=1)))
2106    n=size(a,dim=1)
2107    @:ASSERT(n>0)
2108    if (present(itype)) then
2109      iitype = itype
2110    else
2111      iitype = 1
2112    end if
2113    @:ASSERT(iitype >= 1 .and. iitype <= 3 )
2114
2115    ! Workspace query
2116    allocate(work(1))
2117    allocate(iwork(1))
2118  #:if VTYPE == 'complex'
2119    allocate(rwork(1))
2120  #:endif
2121
2122    call magmaf_${NAME}$_m(ngpus, iitype, jobz, uplo, n, a, n, b, n, w, work, -1,&
2123      #:if VTYPE == 'complex'
2124        & rwork, -1,&
2125      #:endif
2126        & iwork, -1, info)
2127
2128    if (info /= 0) then
2129      call error("Failure in MAGMA_${NAME}$ to determine optimum workspace")
2130    endif
2131
2132 #:if VTYPE == 'complex'
2133    lwork = floor(real(work(1)))
2134    lrwork = floor(rwork(1))
2135    liwork = int(iwork(1))
2136    deallocate(work) ;  allocate(work(lwork))
2137    deallocate(rwork) ;  allocate(rwork(lrwork))
2138    deallocate(iwork) ; allocate(iwork(liwork))
2139  #:endif
2140  #:if VTYPE == 'real'
2141    lwork = floor(work(1))
2142    liwork = int(iwork(1))
2143    deallocate(work) ;  allocate(work(lwork))
2144    deallocate(iwork) ; allocate(iwork(liwork))
2145   #:endif
2146
2147    ! MAGMA Diagonalization
2148    call magmaf_${NAME}$_m(ngpus, iitype, jobz, uplo, n, a, n, b, n, w, work, lwork,&
2149      #:if VTYPE == 'complex'
2150        & rwork, lrwork,&
2151      #:endif
2152        & iwork, liwork, info)
2153
2154    ! test for errors
2155    if (info /= 0) then
2156      if (info < 0) then
2157        write(error_string, "('Failure in diagonalisation routine magmaf_${NAME}$_m,&
2158            & illegal argument at position ',i6)") info
2159        call error(error_string)
2160      else if (info <= n) then
2161        write(error_string, "('Failure in diagonalisation routine magmaf_${NAME}$_m,&
2162            & diagonal element ',i6,' did not converge to zero.')") info
2163        call error(error_string)
2164      else
2165        write(error_string, "('Failure in diagonalisation routine magmaf_${NAME}$_m,&
2166            & non-positive definite overlap! ',i6)") info - n
2167        call error(error_string)
2168      endif
2169    endif
2170
2171  end subroutine ${DTYPE}$_magma_${NAME}$
2172
2173#:endfor
2174
2175#:endif
2176
2177
2178end module dftbp_eigensolver
2179