From dcc9912c989050c99dca96ce535148cfeb42874b Mon Sep 17 00:00:00 2001 From: Rocco Meli Date: Fri, 6 Dec 2024 09:24:05 +0100 Subject: [PATCH] Partial eigenspectrum (#26) --- ci/ci-common.yml | 4 +- ci/docker/build.Dockerfile | 4 +- ci/docker/common.yaml | 3 + spack/packages/dla-future-fortran/package.py | 5 + src/dlaf_fortran.f90 | 461 ++++++++++++++++++- test/helpers/pxheevd.fypp | 81 +++- test/helpers/pxhegvd.fypp | 122 +++-- 7 files changed, 636 insertions(+), 44 deletions(-) diff --git a/ci/ci-common.yml b/ci/ci-common.yml index 0a95e9d..f37be98 100644 --- a/ci/ci-common.yml +++ b/ci/ci-common.yml @@ -30,8 +30,8 @@ stages: reports: dotenv: build.env variables: - SPACK_SHA: 05c7ff4595a4e574cbcf45475c626b9b94c22af1 - SPACK_BUILDCACHE: develop-2024-06-02 + SPACK_SHA: develop-2024-12-01 + SPACK_BUILDCACHE: develop-2024-12-01 SPACK_DLAF_FORTRAN_REPO: ./spack DOCKER_BUILD_ARGS: '[ "BASE_IMAGE", diff --git a/ci/docker/build.Dockerfile b/ci/docker/build.Dockerfile index e377e0b..4aed3f3 100644 --- a/ci/docker/build.Dockerfile +++ b/ci/docker/build.Dockerfile @@ -59,7 +59,9 @@ RUN spack external find \ # Enable Spack build cache ARG SPACK_BUILDCACHE RUN spack mirror add ${SPACK_BUILDCACHE} https://binaries.spack.io/${SPACK_BUILDCACHE} -RUN spack buildcache keys --install --trust --force +RUN spack mirror add develop https://binaries.spack.io/develop && \ + spack buildcache keys --install --trust --force && \ + spack mirror rm develop # Add custom Spack repo ARG SPACK_DLAF_FORTRAN_REPO diff --git a/ci/docker/common.yaml b/ci/docker/common.yaml index 28b8e83..c77cd2e 100644 --- a/ci/docker/common.yaml +++ b/ci/docker/common.yaml @@ -39,6 +39,9 @@ packages: hwloc: variants: - '~libxml2' + dla-future: + require: + - '@master' git: # Force git as non-buildable to allow deprecated versions in environments # https://github.com/spack/spack/pull/30040 diff --git a/spack/packages/dla-future-fortran/package.py b/spack/packages/dla-future-fortran/package.py index 8a92614..608b5a9 100644 --- a/spack/packages/dla-future-fortran/package.py +++ b/spack/packages/dla-future-fortran/package.py @@ -24,6 +24,10 @@ class DlaFutureFortran(CMakePackage): version("0.2.0", sha256="7fd3e1779c111b35f0d2701a024398b4f6e8dea4af523b6c8617d28c0b7ae61a") version("0.1.0", sha256="9fd8a105cbb2f3e1daf8a49910f98fce68ca0b954773dba98a91464cf2e7c1da") + depends_on("c", type="build") + depends_on("cxx", type="build") + depends_on("fortran", type="build") + variant("shared", default=True, description="Build shared libraries.") variant("test", default=False, description="Build tests.") @@ -33,6 +37,7 @@ class DlaFutureFortran(CMakePackage): depends_on("dla-future@0.4.1:0.5 +scalapack", when="@0.1.0") depends_on("dla-future@0.6.0: +scalapack", when="@0.2.0:") + depends_on("dla-future@master: +scalapack", when="@main") depends_on("dla-future +shared", when="+shared") depends_on("mpi", when="+test") diff --git a/src/dlaf_fortran.f90 b/src/dlaf_fortran.f90 index c5fc5bb..f245bc4 100644 --- a/src/dlaf_fortran.f90 +++ b/src/dlaf_fortran.f90 @@ -10,7 +10,7 @@ module dlaf_fortran - use iso_fortran_env, only: dp => real64, sp => real32 + use iso_fortran_env, only: dp => real64, sp => real32, i8 => int64 use iso_c_binding, only: & c_char, & @@ -18,6 +18,7 @@ module dlaf_fortran c_int, & c_loc, & c_ptr, & + c_ptrdiff_t, & c_signed_char, & c_null_char @@ -29,8 +30,14 @@ module dlaf_fortran public :: dlaf_create_grid_from_blacs, dlaf_free_grid public :: dlaf_pspotrf, dlaf_pdpotrf, dlaf_pcpotrf, dlaf_pzpotrf public :: dlaf_pssyevd, dlaf_pdsyevd, dlaf_pcheevd, dlaf_pzheevd + public :: dlaf_pssyevd_partial_spectrum, dlaf_pdsyevd_partial_spectrum + public :: dlaf_pcheevd_partial_spectrum, dlaf_pzheevd_partial_spectrum public :: dlaf_pssygvd, dlaf_pdsygvd, dlaf_pchegvd, dlaf_pzhegvd + public :: dlaf_pssygvd_partial_spectrum, dlaf_pdsygvd_partial_spectrum + public :: dlaf_pchegvd_partial_spectrum, dlaf_pzhegvd_partial_spectrum public :: dlaf_pssygvd_factorized, dlaf_pdsygvd_factorized, dlaf_pchegvd_factorized, dlaf_pzhegvd_factorized + public :: dlaf_pssygvd_partial_spectrum_factorized, dlaf_pdsygvd_partial_spectrum_factorized + public :: dlaf_pchegvd_partial_spectrum_factorized, dlaf_pzhegvd_partial_spectrum_factorized contains @@ -239,6 +246,42 @@ end subroutine dlaf_pssyevd_c end subroutine dlaf_pssyevd + subroutine dlaf_pssyevd_partial_spectrum(uplo, n, a, ia, ja, desca, w, z, iz, jz, descz, il, iu, info) + character, intent(in) :: uplo + integer, intent(in) :: n, ia, ja, iz, jz + integer(kind=i8), intent(in) :: il, iu + integer, dimension(9), intent(in) :: desca, descz + integer, target, intent(out) :: info + real(kind=sp), dimension(:, :), target, intent(inout) :: a, z + real(kind=sp), dimension(:), target, intent(out) :: w + + interface + subroutine dlaf_pssyevd_ps_c(uplo_, n_, a_, ia_, ja_, desca_, w_, z_, iz_, jz_, descz_, il_, iu_, info_) & + bind(C, name='dlaf_pssyevd_partial_spectrum') + + import :: c_int, c_ptr, c_signed_char, c_ptrdiff_t + + integer(kind=c_signed_char), value :: uplo_ + integer(kind=c_int), value :: n_, ia_, ja_, iz_, jz_ + integer(kind=c_ptrdiff_t), value :: il_, iu_ + type(c_ptr), value :: a_, w_, z_ + integer(kind=c_int), dimension(9) :: desca_, descz_ + type(c_ptr), value :: info_ + end subroutine dlaf_pssyevd_ps_c + end interface + + info = -1 + + call dlaf_pssyevd_ps_c(iachar(uplo, c_signed_char), n, & + c_loc(a(1, 1)), ia, ja, desca, & + c_loc(w(1)), & + c_loc(z(1, 1)), iz, jz, descz, & + il - 1, iu, & + c_loc(info) & + ) + + end subroutine dlaf_pssyevd_partial_spectrum + subroutine dlaf_pdsyevd(uplo, n, a, ia, ja, desca, w, z, iz, jz, descz, info) character, intent(in) :: uplo integer, intent(in) :: n, ia, ja, iz, jz @@ -272,6 +315,42 @@ end subroutine dlaf_pdsyevd_c end subroutine dlaf_pdsyevd + subroutine dlaf_pdsyevd_partial_spectrum(uplo, n, a, ia, ja, desca, w, z, iz, jz, descz, il, iu, info) + character, intent(in) :: uplo + integer, intent(in) :: n, ia, ja, iz, jz + integer(kind=i8), intent(in) :: il, iu + integer, dimension(9), intent(in) :: desca, descz + integer, target, intent(out) :: info + real(kind=dp), dimension(:, :), target, intent(inout) :: a, z + real(kind=dp), dimension(:), target, intent(out) :: w + + interface + subroutine dlaf_pdsyevd_ps_c(uplo_, n_, a_, ia_, ja_, desca_, w_, z_, iz_, jz_, descz_, il_, iu_, info_) & + bind(C, name='dlaf_pdsyevd_partial_spectrum') + + import :: c_int, c_ptr, c_signed_char, c_ptrdiff_t + + integer(kind=c_signed_char), value :: uplo_ + integer(kind=c_int), value :: n_, ia_, ja_, iz_, jz_ + integer(kind=c_ptrdiff_t), value :: il_, iu_ + type(c_ptr), value :: a_, w_, z_ + integer(kind=c_int), dimension(9) :: desca_, descz_ + type(c_ptr), value :: info_ + end subroutine dlaf_pdsyevd_ps_c + end interface + + info = -1 + + call dlaf_pdsyevd_ps_c(iachar(uplo, c_signed_char), n, & + c_loc(a(1, 1)), ia, ja, desca, & + c_loc(w(1)), & + c_loc(z(1, 1)), iz, jz, descz, & + il - 1, iu, & + c_loc(info) & + ) + + end subroutine dlaf_pdsyevd_partial_spectrum + subroutine dlaf_pcheevd(uplo, n, a, ia, ja, desca, w, z, iz, jz, descz, info) character, intent(in) :: uplo integer, intent(in) :: n, ia, ja, iz, jz @@ -305,6 +384,42 @@ end subroutine dlaf_pcheevd_c end subroutine dlaf_pcheevd + subroutine dlaf_pcheevd_partial_spectrum(uplo, n, a, ia, ja, desca, w, z, iz, jz, descz, il, iu, info) + character, intent(in) :: uplo + integer, intent(in) :: n, ia, ja, iz, jz + integer(kind=i8), intent(in) :: il, iu + integer, dimension(9), intent(in) :: desca, descz + integer, target, intent(out) :: info + complex(kind=sp), dimension(:, :), target, intent(inout) :: a, z + real(kind=sp), dimension(:), target, intent(out) :: w + + interface + subroutine dlaf_pcheevd_ps_c(uplo_, n_, a_, ia_, ja_, desca_, w_, z_, iz_, jz_, descz_, il_, iu_, info_) & + bind(C, name='dlaf_pcheevd_partial_spectrum') + + import :: c_int, c_ptr, c_signed_char, c_ptrdiff_t + + integer(kind=c_signed_char), value :: uplo_ + integer(kind=c_int), value :: n_, ia_, ja_, iz_, jz_ + integer(kind=c_ptrdiff_t), value :: il_, iu_ + type(c_ptr), value :: a_, w_, z_ + integer(kind=c_int), dimension(9) :: desca_, descz_ + type(c_ptr), value :: info_ + end subroutine dlaf_pcheevd_ps_c + end interface + + info = -1 + + call dlaf_pcheevd_ps_c(iachar(uplo, c_signed_char), n, & + c_loc(a(1, 1)), ia, ja, desca, & + c_loc(w(1)), & + c_loc(z(1, 1)), iz, jz, descz, & + il - 1, iu, & + c_loc(info) & + ) + + end subroutine dlaf_pcheevd_partial_spectrum + subroutine dlaf_pzheevd(uplo, n, a, ia, ja, desca, w, z, iz, jz, descz, info) character, intent(in) :: uplo integer, intent(in) :: n, ia, ja, iz, jz @@ -338,6 +453,42 @@ end subroutine dlaf_pzheevd_c end subroutine dlaf_pzheevd + subroutine dlaf_pzheevd_partial_spectrum(uplo, n, a, ia, ja, desca, w, z, iz, jz, descz, il, iu, info) + character, intent(in) :: uplo + integer, intent(in) :: n, ia, ja, iz, jz + integer(kind=i8), intent(in) :: il, iu + integer, dimension(9), intent(in) :: desca, descz + integer, target, intent(out) :: info + complex(kind=dp), dimension(:, :), target, intent(inout) :: a, z + real(kind=dp), dimension(:), target, intent(out) :: w + + interface + subroutine dlaf_pzheevd_ps_c(uplo_, n_, a_, ia_, ja_, desca_, w_, z_, iz_, jz_, descz_, il_, iu_, info_) & + bind(C, name='dlaf_pzheevd_partial_spectrum') + + import :: c_int, c_ptr, c_signed_char, c_ptrdiff_t + + integer(kind=c_signed_char), value :: uplo_ + integer(kind=c_int), value :: n_, ia_, ja_, iz_, jz_ + integer(kind=c_ptrdiff_t), value :: il_, iu_ + type(c_ptr), value :: a_, w_, z_ + integer(kind=c_int), dimension(9) :: desca_, descz_ + type(c_ptr), value :: info_ + end subroutine dlaf_pzheevd_ps_c + end interface + + info = -1 + + call dlaf_pzheevd_ps_c(iachar(uplo, c_signed_char), n, & + c_loc(a(1, 1)), ia, ja, desca, & + c_loc(w(1)), & + c_loc(z(1, 1)), iz, jz, descz, & + il - 1, iu, & + c_loc(info) & + ) + + end subroutine dlaf_pzheevd_partial_spectrum + subroutine dlaf_pssygvd(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 @@ -372,6 +523,43 @@ end subroutine dlaf_pssygvd_c end subroutine dlaf_pssygvd + subroutine dlaf_pssygvd_partial_spectrum(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, il, iu, info) + character, intent(in) :: uplo + integer, intent(in) :: n, ia, ja, ib, jb, iz, jz + integer(kind=i8), intent(in) :: il, iu + 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_ps_c(uplo_, n_, a_, ia_, ja_, desca_, b_, ib_, jb_, descb_, w_, z_, iz_, jz_, descz_, il_, iu_, info_) & + bind(C, name='dlaf_pssygvd_partial_spectrum') + + import :: c_int, c_ptr, c_signed_char, c_ptrdiff_t + + integer(kind=c_signed_char), value :: uplo_ + integer(kind=c_int), value :: n_, ia_, ja_, ib_, jb_, iz_, jz_ + integer(kind=c_ptrdiff_t), value :: il_, iu_ + 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_ps_c + end interface + + info = -1 + + call dlaf_pssygvd_ps_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, & + il - 1, iu, & + c_loc(info) & + ) + + end subroutine dlaf_pssygvd_partial_spectrum + 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 @@ -408,6 +596,45 @@ end subroutine dlaf_pssygvd_factorized_c end subroutine dlaf_pssygvd_factorized + subroutine dlaf_pssygvd_partial_spectrum_factorized(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, il, iu, info) + character, intent(in) :: uplo + integer, intent(in) :: n, ia, ja, ib, jb, iz, jz + integer(kind=i8), intent(in) :: il, iu + 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_ps_factorized_c( & + uplo_, n_, a_, ia_, ja_, desca_, b_, ib_, jb_, descb_, w_, z_, iz_, jz_, descz_, il_, iu_, info_ & + ) & + bind(C, name='dlaf_pssygvd_partial_spectrum_factorized') + + import :: c_int, c_ptr, c_signed_char, c_ptrdiff_t + + integer(kind=c_signed_char), value :: uplo_ + integer(kind=c_int), value :: n_, ia_, ja_, ib_, jb_, iz_, jz_ + integer(kind=c_ptrdiff_t), value :: il_, iu_ + 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_ps_factorized_c + end interface + + info = -1 + + call dlaf_pssygvd_ps_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, & + il - 1, iu, & + c_loc(info) & + ) + + end subroutine dlaf_pssygvd_partial_spectrum_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 @@ -442,6 +669,43 @@ end subroutine dlaf_pdsygvd_c end subroutine dlaf_pdsygvd + subroutine dlaf_pdsygvd_partial_spectrum(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, il, iu, info) + character, intent(in) :: uplo + integer, intent(in) :: n, ia, ja, ib, jb, iz, jz + integer(kind=i8), intent(in) :: il, iu + 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_ps_c(uplo_, n_, a_, ia_, ja_, desca_, b_, ib_, jb_, descb_, w_, z_, iz_, jz_, descz_, il_, iu_, info_) & + bind(C, name='dlaf_pdsygvd_partial_spectrum') + + import :: c_int, c_ptr, c_signed_char, c_ptrdiff_t + + integer(kind=c_signed_char), value :: uplo_ + integer(kind=c_int), value :: n_, ia_, ja_, ib_, jb_, iz_, jz_ + integer(kind=c_ptrdiff_t), value :: il_, iu_ + 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_ps_c + end interface + + info = -1 + + call dlaf_pdsygvd_ps_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, & + il - 1, iu, & + c_loc(info) & + ) + + end subroutine dlaf_pdsygvd_partial_spectrum + 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 @@ -478,6 +742,45 @@ end subroutine dlaf_pdsygvd_factorized_c end subroutine dlaf_pdsygvd_factorized + subroutine dlaf_pdsygvd_partial_spectrum_factorized(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, il, iu, info) + character, intent(in) :: uplo + integer, intent(in) :: n, ia, ja, ib, jb, iz, jz + integer(kind=i8), intent(in) :: il, iu + 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_ps_factorized_c( & + uplo_, n_, a_, ia_, ja_, desca_, b_, ib_, jb_, descb_, w_, z_, iz_, jz_, descz_, il_, iu_, info_ & + ) & + bind(C, name='dlaf_pdsygvd_partial_spectrum_factorized') + + import :: c_int, c_ptr, c_signed_char, c_ptrdiff_t + + integer(kind=c_signed_char), value :: uplo_ + integer(kind=c_int), value :: n_, ia_, ja_, ib_, jb_, iz_, jz_ + integer(kind=c_ptrdiff_t), value :: il_, iu_ + 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_ps_factorized_c + end interface + + info = -1 + + call dlaf_pdsygvd_ps_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, & + il - 1, iu, & + c_loc(info) & + ) + + end subroutine dlaf_pdsygvd_partial_spectrum_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 @@ -512,6 +815,43 @@ end subroutine dlaf_pchegvd_c end subroutine dlaf_pchegvd + subroutine dlaf_pchegvd_partial_spectrum(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, il, iu, info) + character, intent(in) :: uplo + integer, intent(in) :: n, ia, ja, ib, jb, iz, jz + integer(kind=i8), intent(in) :: il, iu + 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_ps_c(uplo_, n_, a_, ia_, ja_, desca_, b_, ib_, jb_, descb_, w_, z_, iz_, jz_, descz_, il_, iu_, info_) & + bind(C, name='dlaf_pchegvd_partial_spectrum') + + import :: c_int, c_ptr, c_signed_char, c_ptrdiff_t + + integer(kind=c_signed_char), value :: uplo_ + integer(kind=c_int), value :: n_, ia_, ja_, ib_, jb_, iz_, jz_ + integer(kind=c_ptrdiff_t), value :: il_, iu_ + 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_ps_c + end interface + + info = -1 + + call dlaf_pchegvd_ps_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, & + il - 1, iu, & + c_loc(info) & + ) + + end subroutine dlaf_pchegvd_partial_spectrum + 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 @@ -548,6 +888,45 @@ end subroutine dlaf_pchegvd_factorized_c end subroutine dlaf_pchegvd_factorized + subroutine dlaf_pchegvd_partial_spectrum_factorized(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, il, iu, info) + character, intent(in) :: uplo + integer, intent(in) :: n, ia, ja, ib, jb, iz, jz + integer(kind=i8), intent(in) :: il, iu + 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_ps_factorized_c( & + uplo_, n_, a_, ia_, ja_, desca_, b_, ib_, jb_, descb_, w_, z_, iz_, jz_, descz_, il_, iu_, info_ & + ) & + bind(C, name='dlaf_pchegvd_partial_spectrum_factorized') + + import :: c_int, c_ptr, c_signed_char, c_ptrdiff_t + + integer(kind=c_signed_char), value :: uplo_ + integer(kind=c_int), value :: n_, ia_, ja_, ib_, jb_, iz_, jz_ + integer(kind=c_ptrdiff_t), value :: il_, iu_ + 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_ps_factorized_c + end interface + + info = -1 + + call dlaf_pchegvd_ps_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, & + il - 1, iu, & + c_loc(info) & + ) + + end subroutine dlaf_pchegvd_partial_spectrum_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 @@ -582,6 +961,45 @@ end subroutine dlaf_pzhegvd_c end subroutine dlaf_pzhegvd + subroutine dlaf_pzhegvd_partial_spectrum(uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, il, iu, info) + character, intent(in) :: uplo + integer, intent(in) :: n, ia, ja, ib, jb, iz, jz + integer(kind=i8), intent(in) :: il, iu + 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_ps_c( & + uplo_, n_, a_, ia_, ja_, desca_, b_, ib_, jb_, descb_, w_, z_, iz_, jz_, descz_, il_, iu_, info_ & + ) & + bind(C, name='dlaf_pzhegvd_partial_spectrum') + + import :: c_int, c_ptr, c_signed_char, c_ptrdiff_t + + integer(kind=c_signed_char), value :: uplo_ + integer(kind=c_int), value :: n_, ia_, ja_, ib_, jb_, iz_, jz_ + integer(kind=c_ptrdiff_t), value :: il_, iu_ + 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_ps_c + end interface + + info = -1 + + call dlaf_pzhegvd_ps_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, & + il - 1, iu, & + c_loc(info) & + ) + + end subroutine dlaf_pzhegvd_partial_spectrum + 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 @@ -618,4 +1036,45 @@ end subroutine dlaf_pzhegvd_factorized_c end subroutine dlaf_pzhegvd_factorized + subroutine dlaf_pzhegvd_partial_spectrum_factorized( & + uplo, n, a, ia, ja, desca, b, ib, jb, descb, w, z, iz, jz, descz, il, iu, info & + ) + character, intent(in) :: uplo + integer, intent(in) :: n, ia, ja, ib, jb, iz, jz + integer(kind=i8), intent(in) :: il, iu + 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_ps_factorized_c( & + uplo_, n_, a_, ia_, ja_, desca_, b_, ib_, jb_, descb_, w_, z_, iz_, jz_, descz_, il_, iu_, info_ & + ) & + bind(C, name='dlaf_pzhegvd_partial_spectrum_factorized') + + import :: c_int, c_ptr, c_signed_char, c_ptrdiff_t + + integer(kind=c_signed_char), value :: uplo_ + integer(kind=c_int), value :: n_, ia_, ja_, ib_, jb_, iz_, jz_ + integer(kind=c_ptrdiff_t), value :: il_, iu_ + 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_ps_factorized_c + end interface + + info = -1 + + call dlaf_pzhegvd_ps_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, & + il - 1, iu, & + c_loc(info) & + ) + + end subroutine dlaf_pzhegvd_partial_spectrum_factorized + end module dlaf_fortran diff --git a/test/helpers/pxheevd.fypp b/test/helpers/pxheevd.fypp index 8a65fb3..5c50501 100644 --- a/test/helpers/pxheevd.fypp +++ b/test/helpers/pxheevd.fypp @@ -13,12 +13,12 @@ #:set names = {('sp', 'real'): 'ssy', ('sp', 'complex'): 'che', ('dp', 'real'): 'dsy', ('dp', 'complex'): 'zhe'} #:set symbols = {('sp', 'real'): 's', ('sp', 'complex'): 'c', ('dp', 'real'): 'd', ('dp', 'complex'): 'z'} module pxheevd_tests - use iso_fortran_env, only: error_unit, sp => real32, dp => real64 + use iso_fortran_env, only: error_unit, sp => real32, dp => real64, i8 => int64 use dlaf_fortran, only: dlaf_initialize, dlaf_finalize, dlaf_create_grid_from_blacs, dlaf_free_grid #:for dtype in precision #:for type in types #:set name = names[(dtype, type)] - use dlaf_fortran, only: dlaf_p${name}$evd + use dlaf_fortran, only: dlaf_p${name}$evd, dlaf_p${name}$evd_partial_spectrum #:endfor #:endfor @@ -50,6 +50,8 @@ contains subroutine p${name}$evd_test integer, parameter :: n = 4 + integer(kind=i8), parameter :: eval_idx_begin = 1 + integer(kind=i8), parameter :: eval_idx_end = 2 integer:: nprow, npcol integer:: i, j @@ -59,12 +61,15 @@ contains integer :: ictxt, ictxt_0 integer :: info, lld, nb, ma, na integer :: desca(9), desca_local_dlaf(9), desca_local_scalapack(9) - integer :: descz_local_dlaf(9), descz_local_scalapack(9) - integer :: descz_dlaf(9), descz_scalapack(9) + integer :: descz_local_dlaf(9), descz_local_scalapack(9), descz_local_dlaf_partial_spectrum(9) + integer :: descz_dlaf(9), descz_scalapack(9), descz_dlaf_partial_spectrum(9) ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: A, A_local_dlaf, A_local_scalapack + ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: A_local_dlaf_partial_spectrum ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: Z_local_dlaf, Z_local_scalapack + ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: Z_local_dlaf_partial_spectrum ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: Z_dlaf, Z_scalapack - real(kind=${dtype}$), dimension(:), allocatable :: W_dlaf, W_scalapack + ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: Z_dlaf_partial_spectrum + real(kind=${dtype}$), dimension(:), allocatable :: W_dlaf, W_scalapack, W_dlaf_partial_spectrum real(kind=${dtype}$), parameter :: abstol = #{if dtype == 'sp'}# 1e-5_${dtype}$ #{else}# 1e-8_${dtype}$ #{endif}# integer, parameter :: lwork = 100, liwork = 100 @@ -93,14 +98,17 @@ contains call init_desc(desca) call init_desc(descz_scalapack) call init_desc(descz_dlaf) + call init_desc(descz_dlaf_partial_spectrum) if (rank == 0) then allocate (A(n, n)) allocate (Z_dlaf(n, n)) allocate (Z_scalapack(n, n)) + allocate (Z_dlaf_partial_spectrum(n, n)) call descinit(desca, n, n, n, n, 0, 0, ictxt_0, n, info) call descinit(descz_dlaf, n, n, n, n, 0, 0, ictxt_0, n, info) call descinit(descz_scalapack, n, n, n, n, 0, 0, ictxt_0, n, info) + call descinit(descz_dlaf_partial_spectrum, n, n, n, n, 0, 0, ictxt_0, n, info) call set_random_matrix(A) end if @@ -109,17 +117,20 @@ 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 (Z_local_dlaf(ma, na), Z_local_scalapack(ma, na)) - allocate (W_dlaf(n), W_scalapack(n)) + allocate (A_local_dlaf(ma, na), A_local_scalapack(ma, na), A_local_dlaf_partial_spectrum(ma, na)) + allocate (Z_local_dlaf(ma, na), Z_local_scalapack(ma, na), Z_local_dlaf_partial_spectrum(ma, na)) + allocate (W_dlaf(n), W_scalapack(n), W_dlaf_partial_spectrum(n)) ! + ---- + ! | DLAF | ! + ---- + call descinit(desca_local_dlaf, n, n, nb, nb, 0, 0, ictxt, lld, info) - call descinit(descz_local_dlaf, n, n, nb, nb, 0, 0, ictxt, lld, info) 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, a, 1, 1, desca, A_local_dlaf_partial_spectrum, 1, 1, desca_local_dlaf, ictxt) + + call descinit(descz_local_dlaf, n, n, nb, nb, 0, 0, ictxt, lld, info) + call descinit(descz_local_dlaf_partial_spectrum, n, n, nb, nb, 0, 0, ictxt, lld, info) ! Solve with DLAF call dlaf_initialize() @@ -130,6 +141,13 @@ contains W_dlaf, Z_local_dlaf, 1, 1, descz_local_dlaf, & info & ) + call dlaf_p${name}$evd_partial_spectrum( & + 'L', & + n, A_local_dlaf_partial_spectrum, 1, 1, desca_local_dlaf, & + W_dlaf_partial_spectrum, Z_local_dlaf_partial_spectrum, 1, 1, descz_local_dlaf_partial_spectrum, & + eval_idx_begin, eval_idx_end, & + info & + ) call dlaf_free_grid(ictxt) call dlaf_finalize() if (info /= 0) then @@ -138,6 +156,12 @@ 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_partial_spectrum, 1, 1, descz_local_dlaf_partial_spectrum, & + Z_dlaf_partial_spectrum, 1, 1, descz_dlaf_partial_spectrum, & + ictxt & + ) ! + --------- + ! | ScaLAPACK | @@ -170,7 +194,9 @@ contains call p${symbol}$gemr2d(n, n, Z_local_scalapack, 1, 1, descz_local_scalapack, Z_scalapack, 1, 1, descz_scalapack, ictxt) - ! Check results + ! + ------------- + + ! | Check Results | + ! + ------------- + ! Results are checked only on rank 0 failed = .false. @@ -179,6 +205,11 @@ contains failed = .true. write (error_unit, *) "ERROR: DLAF != ScaLAPACK (eigenvalues)" end if + + if (.not. allclose(W_dlaf_partial_spectrum, W_scalapack)) then + failed = .true. + write (error_unit, *) "ERROR: DLAF_PS != ScaLAPACK (eigenvalues)" + end if end if call bcast_check(failed) @@ -188,10 +219,32 @@ contains failed = .false. if (rank == 0) then - if (.not. allclose(abs(Z_dlaf), abs(Z_scalapack), atol=abstol)) then - failed = .true. - write (error_unit, *) "ERROR: DLAF != ScaLAPACK (eigenvectors)" - end if + if (.not. allclose(abs(Z_dlaf), abs(Z_scalapack), atol=abstol)) then + failed = .true. + write (error_unit, *) "ERROR: DLAF != ScaLAPACK (eigenvectors)" + end if + + ! The first "eval_idx_end" eigenvectors are the same as ScaLAPACK + if ( & + .not. allclose( & + abs(Z_dlaf_partial_spectrum(:,eval_idx_begin:eval_idx_end)), & + abs(Z_scalapack(:, eval_idx_begin:eval_idx_end)), & + atol=abstol) & + ) then + failed = .true. + write (error_unit, *) "ERROR: DLAF_PS != ScaLAPACK (eigenvectors [begin,end])" + end if + + ! Something is off if the lastt eigenvectors are the same as ScaLAPACK + if ( & + allclose( & + abs(Z_dlaf_partial_spectrum(:,eval_idx_end+1:n)), & + abs(Z_scalapack(:,eval_idx_end+1:n)), & + atol=abstol) & + ) then + failed = .true. + write (error_unit, *) "ERROR: DLAF_PS == ScaLAPACK (eigenvectors [end+1:N])" + end if end if call bcast_check(failed) diff --git a/test/helpers/pxhegvd.fypp b/test/helpers/pxhegvd.fypp index 4ac02ae..b82147f 100644 --- a/test/helpers/pxhegvd.fypp +++ b/test/helpers/pxhegvd.fypp @@ -13,13 +13,15 @@ #:set names = {('sp', 'real'): 'ssy', ('sp', 'complex'): 'che', ('dp', 'real'): 'dsy', ('dp', 'complex'): 'zhe'} #:set symbols = {('sp', 'real'): 's', ('sp', 'complex'): 'c', ('dp', 'real'): 'd', ('dp', 'complex'): 'z'} module pxhegvd_tests - use iso_fortran_env, only: error_unit, sp => real32, dp => real64 + use iso_fortran_env, only: error_unit, sp => real32, dp => real64, i8 => int64 use dlaf_fortran, only: dlaf_initialize, dlaf_finalize, dlaf_create_grid_from_blacs, dlaf_free_grid #:for dtype in precision #: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_partial_spectrum use dlaf_fortran, only: dlaf_p${name}$gvd_factorized + use dlaf_fortran, only: dlaf_p${name}$gvd_partial_spectrum_factorized #:endfor #:endfor @@ -53,6 +55,8 @@ contains subroutine p${name}$gvd_test integer, parameter :: n = 4 + integer(kind=i8), parameter :: eval_idx_begin = 1 + integer(kind=i8), parameter :: eval_idx_end = 2 integer:: nprow, npcol integer:: i, j @@ -70,10 +74,13 @@ contains 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, A_local_store - ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: B, B_local_dlaf, B_local_scalapack + ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: B, B_local_dlaf, B_local_scalapack, B_local_store ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: Z_local_dlaf, Z_local_scalapack, Z_local_dlaf_factorized + ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: Z_local_dlaf_partial_spectrum, Z_local_dlaf_partial_spectrum_factorized ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: Z_dlaf, Z_scalapack, Z_dlaf_factorized + ${type}$ (kind=${dtype}$), dimension(:, :), allocatable :: Z_dlaf_partial_spectrum, Z_dlaf_partial_spectrum_factorized real(kind=${dtype}$), dimension(:), allocatable :: W_dlaf, W_scalapack, W_dlaf_factorized + real(kind=${dtype}$), dimension(:), allocatable :: W_dlaf_partial_spectrum, W_dlaf_partial_spectrum_factorized real(kind=${dtype}$), parameter :: abstol = #{if dtype == 'sp'}# 1e-5_${dtype}$ #{else}# 1e-8_${dtype}$ #{endif}# integer, parameter :: lwork = 100, liwork = 100 @@ -112,6 +119,8 @@ contains allocate (Z_dlaf(n, n)) allocate (Z_scalapack(n, n)) allocate (Z_dlaf_factorized(n, n)) + allocate (Z_dlaf_partial_spectrum(n, n)) + allocate (Z_dlaf_partial_spectrum_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) @@ -127,9 +136,11 @@ contains na = numroc(n, nb, mycol, 0, npcol) lld = max(1, ma) 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 (B_local_dlaf(ma, na), B_local_scalapack(ma, na), B_local_store(ma, na)) allocate (Z_local_dlaf(ma, na), Z_local_scalapack(ma, na), Z_local_dlaf_factorized(ma, na)) + allocate (Z_local_dlaf_partial_spectrum(ma, na), Z_local_dlaf_partial_spectrum_factorized(ma, na)) allocate (W_dlaf(n), W_scalapack(n), W_dlaf_factorized(n)) + allocate (W_dlaf_partial_spectrum(n), W_dlaf_partial_spectrum_factorized(n)) ! + ---- + ! | DLAF | @@ -137,12 +148,15 @@ contains call descinit(desca_local_dlaf, n, n, nb, nb, 0, 0, ictxt, lld, info) call descinit(descb_local_dlaf, n, n, nb, nb, 0, 0, ictxt, lld, info) - call descinit(descz_local_dlaf, n, n, nb, nb, 0, 0, ictxt, lld, info) 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) + + call descinit(descz_local_dlaf, n, n, nb, nb, 0, 0, ictxt, lld, info) + call descinit(descz_local_dlaf, n, n, nb, nb, 0, 0, ictxt, lld, info) ! Store local copy A_local_store(:, :) = A_local_dlaf(:, :) + B_local_store(:, :) = B_local_dlaf(:, :) ! Solve with DLAF call dlaf_initialize() @@ -166,6 +180,28 @@ contains W_dlaf_factorized, Z_local_dlaf_factorized, 1, 1, descz_local_dlaf, & info & ) + + ! Solve with Cholesky factorization + A_local_dlaf(:, :) = A_local_store(:, :) + B_local_dlaf(:, :) = B_local_store(:, :) + call dlaf_p${name}$gvd_partial_spectrum( & + 'L', & + n, A_local_dlaf, 1, 1, desca_local_dlaf, & + B_local_dlaf, 1, 1, descb_local_dlaf, & + W_dlaf_partial_spectrum, Z_local_dlaf_partial_spectrum, 1, 1, descz_local_dlaf, & + eval_idx_begin, eval_idx_end, & + 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_partial_spectrum_factorized, Z_local_dlaf_partial_spectrum_factorized, 1, 1, descz_local_dlaf, & + info & + ) call dlaf_free_grid(ictxt) call dlaf_finalize() @@ -176,6 +212,8 @@ contains 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) + call p${symbol}$gemr2d(n, n, Z_local_dlaf_partial_spectrum, 1, 1, descz_local_dlaf, Z_dlaf_partial_spectrum, 1, 1, descz_dlaf, ictxt) + call p${symbol}$gemr2d(n, n, Z_local_dlaf_partial_spectrum_factorized, 1, 1, descz_local_dlaf, Z_dlaf_partial_spectrum_factorized, 1, 1, descz_dlaf, ictxt) ! + --------- + ! | ScaLAPACK | @@ -214,7 +252,9 @@ contains call p${symbol}$gemr2d(n, n, Z_local_scalapack, 1, 1, descz_local_scalapack, Z_scalapack, 1, 1, descz_scalapack, ictxt) - ! Check results + ! + ------------- + + ! | Check Results | + ! + ------------- + ! Results are checked only on rank 0 failed = .false. @@ -223,18 +263,20 @@ contains failed = .true. write (error_unit, *) "ERROR: DLAF != ScaLAPACK (eigenvalues)" 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 + + if (.not. allclose(W_dlaf_factorized, W_scalapack)) then + failed = .true. + write (error_unit, *) "ERROR: DLAF Factorized != ScaLAPACK (eigenvalues)" + end if + + if (.not. allclose(W_dlaf_partial_spectrum, W_scalapack)) then + failed = .true. + write (error_unit, *) "ERROR: DLAF Partial Spectrum != ScaLAPACK (eigenvalues)" + end if + + if (.not. allclose(W_dlaf_partial_spectrum_factorized, W_scalapack)) then failed = .true. - write (error_unit, *) "ERROR: DLAF != DLAF Factorized (eigenvalues)" + write (error_unit, *) "ERROR: DLAF Partial Spectrum Factorized != ScaLAPACK (eigenvalues)" end if end if @@ -249,18 +291,46 @@ contains failed = .true. write (error_unit, *) "ERROR: DLAF != ScaLAPACK (eigenvectors)" 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 + if (.not. allclose(abs(Z_dlaf_factorized), abs(Z_dlaf))) then + failed = .true. + write (error_unit, *) "ERROR: DLAF Factorized != DLAF (eigenvectors)" + end if + + if ( & + .not. allclose( & + abs(Z_dlaf_partial_spectrum(:,eval_idx_begin:eval_idx_end)), & + abs(Z_dlaf(:, eval_idx_begin:eval_idx_end))) & + ) then + failed = .true. + write (error_unit, *) "ERROR: DLAF Partial Spectrum != DLAF (eigenvectors [begin, end])" + end if + + if ( & + allclose( & + abs(Z_dlaf_partial_spectrum(:,eval_idx_end+1:n)), & + abs(Z_dlaf(:,eval_idx_end+1:n))) & + ) then + failed = .true. + write (error_unit, *) "ERROR: DLAF Partial Spectrum != DLAF (eigenvectors [end+1, N])" + end if + + if ( & + .not. allclose( & + abs(Z_dlaf_partial_spectrum_factorized(:,eval_idx_begin:eval_idx_end)), & + abs(Z_dlaf(:, eval_idx_begin:eval_idx_end))) & + ) then + failed = .true. + write (error_unit, *) "ERROR: DLAF Partial Spectrum Factorized != DLAF (eigenvectors [begin, end])" + end if + + if ( & + allclose( & + abs(Z_dlaf_partial_spectrum(:,eval_idx_end+1:n)), & + abs(Z_dlaf(:,eval_idx_end+1:n))) & + ) then failed = .true. - write (error_unit, *) "ERROR: DLAF != DLAF Factorized (eigenvectors)" + write (error_unit, *) "ERROR: DLAF Partial Spectrum Factorized != DLAF (eigenvectors [end+1, N])" end if end if