Generalised eigenvalues solver with pre-factorised B matrix

CHANGELOG.md

@@ -2,6 +2,10 @@
 
 ## DLA-Future-Fortran X.Y.Z
 
+### Added
+
+* Generalized eigensolver for an already factorized $\mathbf{B}$ matrix [PR #18]
+
 ### Changed
 
 * Name of the generalized eigenvalue solver from `*gvx` to `*gvd` [PR #16]

src/dlaf_fortran.f90

@@ -30,6 +30,7 @@ module dlaf_fortran
    public :: dlaf_pspotrf, dlaf_pdpotrf, dlaf_pcpotrf, dlaf_pzpotrf
    public :: dlaf_pssyevd, dlaf_pdsyevd, dlaf_pcheevd, dlaf_pzheevd
    public :: dlaf_pssygvd, dlaf_pdsygvd, dlaf_pchegvd, dlaf_pzhegvd
+   public :: dlaf_pssygvd_factorized, dlaf_pdsygvd_factorized, dlaf_pchegvd_factorized, dlaf_pzhegvd_factorized
 
 contains
 
@@ -371,6 +372,42 @@ end subroutine dlaf_pssygvd_c
 
    end subroutine dlaf_pssygvd
 
+   subroutine dlaf_pssygvd_factorized(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, info)
+      character, intent(in) :: uplo
+      integer, intent(in) :: n, ia, ja, ib, jb, iz, jz
+      integer, dimension(9), intent(in) :: desca, descb, descz
+      integer, target, intent(out) :: info
+      real(kind=sp), dimension(:, :), target, intent(inout) :: a, b, z
+      real(kind=sp), dimension(:), target, intent(out) :: w
+
+      interface
+         subroutine dlaf_pssygvd_factorized_c( &
+            uplo_, n_, a_, ia_, ja_, desca_, b_, ib_, jb_, descb_, w_, z_, iz_, jz_, descz_, info_ &
+            ) &
+            bind(C, name='dlaf_pssygvd_factorized')
+
+            import :: c_int, c_ptr, c_signed_char
+
+            integer(kind=c_signed_char), value :: uplo_
+            integer(kind=c_int), value :: n_, ia_, ja_, ib_, jb_, iz_, jz_
+            type(c_ptr), value :: a_, b_, w_, z_
+            integer(kind=c_int), dimension(9) :: desca_, descb_, descz_
+            type(c_ptr), value :: info_
+         end subroutine dlaf_pssygvd_factorized_c
+      end interface
+
+      info = -1
+
+      call dlaf_pssygvd_factorized_c(iachar(uplo, c_signed_char), n, &
+                                     c_loc(a(1, 1)), ia, ja, desca, &
+                                     c_loc(b(1, 1)), ib, jb, descb, &
+                                     c_loc(w(1)), &
+                                     c_loc(z(1, 1)), iz, jz, descz, &
+                                     c_loc(info) &
+                                     )
+
+   end subroutine dlaf_pssygvd_factorized
+
    subroutine dlaf_pdsygvd(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, info)
       character, intent(in) :: uplo
       integer, intent(in) :: n, ia, ja, ib, jb, iz, jz
@@ -405,6 +442,42 @@ end subroutine dlaf_pdsygvd_c
 
    end subroutine dlaf_pdsygvd
 
+   subroutine dlaf_pdsygvd_factorized(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, info)
+      character, intent(in) :: uplo
+      integer, intent(in) :: n, ia, ja, ib, jb, iz, jz
+      integer, dimension(9), intent(in) :: desca, descb, descz
+      integer, target, intent(out) :: info
+      real(kind=dp), dimension(:, :), target, intent(inout) :: a, b, z
+      real(kind=dp), dimension(:), target, intent(out) :: w
+
+      interface
+         subroutine dlaf_pdsygvd_factorized_c( &
+            uplo_, n_, a_, ia_, ja_, desca_, b_, ib_, jb_, descb_, w_, z_, iz_, jz_, descz_, info_ &
+            ) &
+            bind(C, name='dlaf_pdsygvd_factorized')
+
+            import :: c_int, c_ptr, c_signed_char
+
+            integer(kind=c_signed_char), value :: uplo_
+            integer(kind=c_int), value :: n_, ia_, ja_, ib_, jb_, iz_, jz_
+            type(c_ptr), value :: a_, b_, w_, z_
+            integer(kind=c_int), dimension(9) :: desca_, descb_, descz_
+            type(c_ptr), value :: info_
+         end subroutine dlaf_pdsygvd_factorized_c
+      end interface
+
+      info = -1
+
+      call dlaf_pdsygvd_factorized_c(iachar(uplo, c_signed_char), n, &
+                                     c_loc(a(1, 1)), ia, ja, desca, &
+                                     c_loc(b(1, 1)), ib, jb, descb, &
+                                     c_loc(w(1)), &
+                                     c_loc(z(1, 1)), iz, jz, descz, &
+                                     c_loc(info) &
+                                     )
+
+   end subroutine dlaf_pdsygvd_factorized
+
    subroutine dlaf_pchegvd(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, info)
       character, intent(in) :: uplo
       integer, intent(in) :: n, ia, ja, ib, jb, iz, jz
@@ -439,6 +512,42 @@ end subroutine dlaf_pchegvd_c
 
    end subroutine dlaf_pchegvd
 
+   subroutine dlaf_pchegvd_factorized(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, info)
+      character, intent(in) :: uplo
+      integer, intent(in) :: n, ia, ja, ib, jb, iz, jz
+      integer, dimension(9), intent(in) :: desca, descb, descz
+      integer, target, intent(out) :: info
+      complex(kind=sp), dimension(:, :), target, intent(inout) :: a, b, z
+      real(kind=sp), dimension(:), target, intent(out) :: w
+
+      interface
+         subroutine dlaf_pchegvd_factorized_c( &
+            uplo_, n_, a_, ia_, ja_, desca_, b_, ib_, jb_, descb_, w_, z_, iz_, jz_, descz_, info_ &
+            ) &
+            bind(C, name='dlaf_pchegvd_factorized')
+
+            import :: c_int, c_ptr, c_signed_char
+
+            integer(kind=c_signed_char), value :: uplo_
+            integer(kind=c_int), value :: n_, ia_, ja_, ib_, jb_, iz_, jz_
+            type(c_ptr), value :: a_, b_, w_, z_
+            integer(kind=c_int), dimension(9) :: desca_, descb_, descz_
+            type(c_ptr), value :: info_
+         end subroutine dlaf_pchegvd_factorized_c
+      end interface
+
+      info = -1
+
+      call dlaf_pchegvd_factorized_c(iachar(uplo, c_signed_char), n, &
+                                     c_loc(a(1, 1)), ia, ja, desca, &
+                                     c_loc(b(1, 1)), ib, jb, descb, &
+                                     c_loc(w(1)), &
+                                     c_loc(z(1, 1)), iz, jz, descz, &
+                                     c_loc(info) &
+                                     )
+
+   end subroutine dlaf_pchegvd_factorized
+
    subroutine dlaf_pzhegvd(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, info)
       character, intent(in) :: uplo
       integer, intent(in) :: n, ia, ja, ib, jb, iz, jz
@@ -473,4 +582,40 @@ end subroutine dlaf_pzhegvd_c
 
    end subroutine dlaf_pzhegvd
 
+   subroutine dlaf_pzhegvd_factorized(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, info)
+      character, intent(in) :: uplo
+      integer, intent(in) :: n, ia, ja, ib, jb, iz, jz
+      integer, dimension(9), intent(in) :: desca, descb, descz
+      integer, target, intent(out) :: info
+      complex(kind=dp), dimension(:, :), target, intent(inout) :: a, b, z
+      real(kind=dp), dimension(:), target, intent(out) :: w
+
+      interface
+         subroutine dlaf_pzhegvd_factorized_c( &
+            uplo_, n_, a_, ia_, ja_, desca_, b_, ib_, jb_, descb_, w_, z_, iz_, jz_, descz_, info_ &
+            ) &
+            bind(C, name='dlaf_pzhegvd_factorized')
+
+            import :: c_int, c_ptr, c_signed_char
+
+            integer(kind=c_signed_char), value :: uplo_
+            integer(kind=c_int), value :: n_, ia_, ja_, ib_, jb_, iz_, jz_
+            type(c_ptr), value :: a_, b_, w_, z_
+            integer(kind=c_int), dimension(9) :: desca_, descb_, descz_
+            type(c_ptr), value :: info_
+         end subroutine dlaf_pzhegvd_factorized_c
+      end interface
+
+      info = -1
+
+      call dlaf_pzhegvd_factorized_c(iachar(uplo, c_signed_char), n, &
+                                     c_loc(a(1, 1)), ia, ja, desca, &
+                                     c_loc(b(1, 1)), ib, jb, descb, &
+                                     c_loc(w(1)), &
+                                     c_loc(z(1, 1)), iz, jz, descz, &
+                                     c_loc(info) &
+                                     )
+
+   end subroutine dlaf_pzhegvd_factorized
+
 end module dlaf_fortran

test/helpers/pxhegvd.fypp

@@ -19,6 +19,7 @@ module pxhegvd_tests
       #:for type in types
          #:set name = names[(dtype, type)]
          use dlaf_fortran, only: dlaf_p${name}$gvd
+         use dlaf_fortran, only: dlaf_p${name}$gvd_factorized
       #:endfor
    #:endfor
 
@@ -68,11 +69,11 @@ contains
             integer :: descb(9), descb_local_dlaf(9), descb_local_scalapack(9)
             integer :: descz_local_dlaf(9), descz_local_scalapack(9)
             integer :: descz_dlaf(9), descz_scalapack(9)
-            ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: A, A_local_dlaf, A_local_scalapack
+            ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: A, A_local_dlaf, A_local_scalapack, A_local_store
             ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: B, B_local_dlaf, B_local_scalapack
-            ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: Z_local_dlaf, Z_local_scalapack
-            ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: Z_dlaf, Z_scalapack
-            real(kind=${dtype}$), dimension(:), allocatable :: W_dlaf, W_scalapack
+            ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: Z_local_dlaf, Z_local_scalapack, Z_local_dlaf_factorized
+            ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: Z_dlaf, Z_scalapack, Z_dlaf_factorized
+            real(kind=${dtype}$), dimension(:), allocatable :: W_dlaf, W_scalapack, W_dlaf_factorized
             real(kind=${dtype}$), parameter :: abstol = #{if dtype == 'sp'}# 1e-5_${dtype}$ #{else}# 1e-8_${dtype}$ #{endif}#
 
             integer, parameter :: lwork = 100, liwork = 100
@@ -110,6 +111,7 @@ contains
                allocate (B(n, n))
                allocate (Z_dlaf(n, n))
                allocate (Z_scalapack(n, n))
+               allocate (Z_dlaf_factorized(n, n))
 
                call descinit(desca, n, n, n, n, 0, 0, ictxt_0, n, info)
                call descinit(descb, n, n, n, n, 0, 0, ictxt_0, n, info)
@@ -124,10 +126,10 @@ contains
             ma = numroc(n, nb, myrow, 0, nprow)
             na = numroc(n, nb, mycol, 0, npcol)
             lld = max(1, ma)
-            allocate (A_local_dlaf(ma, na), A_local_scalapack(ma, na))
+            allocate (A_local_dlaf(ma, na), A_local_scalapack(ma, na), A_local_store(ma, na))
             allocate (B_local_dlaf(ma, na), B_local_scalapack(ma, na))
-            allocate (Z_local_dlaf(ma, na), Z_local_scalapack(ma, na))
-            allocate (W_dlaf(n), W_scalapack(n))
+            allocate (Z_local_dlaf(ma, na), Z_local_scalapack(ma, na), Z_local_dlaf_factorized(ma, na))
+            allocate (W_dlaf(n), W_scalapack(n), W_dlaf_factorized(n))
 
             ! + ---- +
             ! | DLAF |
@@ -139,16 +141,32 @@ contains
             call p${symbol}$gemr2d(n, n, a, 1, 1, desca, A_local_dlaf, 1, 1, desca_local_dlaf, ictxt)
             call p${symbol}$gemr2d(n, n, b, 1, 1, desca, B_local_dlaf, 1, 1, descb_local_dlaf, ictxt)
 
+            ! Store local copy
+            A_local_store(:, :) = A_local_dlaf(:, :)
+
             ! Solve with DLAF
             call dlaf_initialize()
             call dlaf_create_grid_from_blacs(ictxt)
+
+            ! Solve with Cholesky factorization
             call dlaf_p${name}$gvd( &
                'L', &
                n, A_local_dlaf, 1, 1, desca_local_dlaf, &
                B_local_dlaf, 1, 1, descb_local_dlaf, &
                W_dlaf, Z_local_dlaf, 1, 1, descz_local_dlaf, &
                info &
                )
+
+            ! Solve with Cholesky factorisation from the previous call 
+            A_local_dlaf(:, :) = A_local_store(:, :)
+            call dlaf_p${name}$gvd_factorized( &
+               'L', &
+               n, A_local_dlaf, 1, 1, desca_local_dlaf, &
+               B_local_dlaf, 1, 1, descb_local_dlaf, &
+               W_dlaf_factorized, Z_local_dlaf_factorized, 1, 1, descz_local_dlaf, &
+               info &
+               )
+
             call dlaf_free_grid(ictxt)
             call dlaf_finalize()
             if (info /= 0) then
@@ -157,6 +175,7 @@ contains
             end if
 
             call p${symbol}$gemr2d(n, n, Z_local_dlaf, 1, 1, descz_local_dlaf, Z_dlaf, 1, 1, descz_dlaf, ictxt)
+            call p${symbol}$gemr2d(n, n, Z_local_dlaf_factorized, 1, 1, descz_local_dlaf, Z_dlaf_factorized, 1, 1, descz_dlaf, ictxt)
 
             ! + --------- +
             ! | ScaLAPACK |
@@ -206,6 +225,19 @@ contains
                end if
             end if
 
+            call bcast_check(failed)
+            if (failed) then
+               call terminate(ictxt)
+            end if
+
+            failed = .false.
+            if (rank == 0) then
+               if (.not. allclose(W_dlaf, W_dlaf_factorized)) then
+                  failed = .true.
+                  write (error_unit, *) "ERROR: DLAF != DLAF Factorized (eigenvalues)"
+               end if
+            end if
+
             call bcast_check(failed)
             if (failed) then
                call terminate(ictxt)
@@ -219,6 +251,19 @@ contains
                end if
             end if
 
+            call bcast_check(failed)
+            if (failed) then
+               call terminate(ictxt)
+            end if
+
+            failed = .false.
+            if (rank == 0) then
+               if (.not. allclose(abs(Z_dlaf), abs(Z_dlaf_factorized))) then
+                  failed = .true.
+                  write (error_unit, *) "ERROR: DLAF != DLAF Factorized (eigenvectors)"
+               end if
+            end if
+
             call bcast_check(failed)
             if (failed) then
                call terminate(ictxt)
@@ -234,12 +279,15 @@ contains
             end if
             if (allocated(A_local_dlaf)) deallocate (A_local_dlaf)
             if (allocated(A_local_scalapack)) deallocate (A_local_scalapack)
+            if (allocated(A_local_store)) deallocate (A_local_store)
             if (allocated(B_local_dlaf)) deallocate (B_local_dlaf)
             if (allocated(B_local_scalapack)) deallocate (B_local_scalapack)
             if (allocated(Z_local_dlaf)) deallocate (Z_local_dlaf)
             if (allocated(Z_local_scalapack)) deallocate (Z_local_scalapack)
+            if (allocated(Z_local_dlaf_factorized)) deallocate (Z_local_dlaf_factorized)
             if (allocated(W_dlaf)) deallocate (W_dlaf)
             if (allocated(W_scalapack)) deallocate (W_scalapack)
+            if (allocated(W_dlaf_factorized)) deallocate (W_dlaf_factorized)
             if (allocated(iclustr)) deallocate (iclustr)
             if (allocated(gap)) deallocate (gap)