generalized eigenvector problem progress

generate_automake_test_programs.py

@@ -133,6 +133,10 @@ for m, g, t, p, d, s, l in product(
     if (t == "qr" and (s == "1stage" or d == "complex")):
         continue
 
+    #TODO: this does not work at the moment
+    if(t == "generalized" and (l == "all_layouts" or s == "2stage")):
+        continue
+
     create_test(m, g, t, p, d, s, l, "fortran")
 
 

src/elpa_impl.F90

@@ -125,12 +125,16 @@ module elpa_impl
      procedure, public :: associate_int => elpa_associate_int  !< public method to set some pointers
 
      procedure, private :: elpa_transform_generalized_d
+     procedure, private :: elpa_transform_back_generalized_d
      procedure, private :: elpa_transform_generalized_dc
+     procedure, private :: elpa_transform_back_generalized_dc
 #ifdef WANT_SINGLE_PRECISION_REAL
      procedure, private :: elpa_transform_generalized_f
+     procedure, private :: elpa_transform_back_generalized_f
 #endif
 #ifdef WANT_SINGLE_PRECISION_COMPLEX
      procedure, private :: elpa_transform_generalized_fc
+     procedure, private :: elpa_transform_back_generalized_fc
 #endif
   end type elpa_impl_t
 
@@ -1236,6 +1240,11 @@ module elpa_impl
       logical             :: success_l
 
       call self%elpa_transform_generalized_d(a, b, sc_desc, error_l)
+      if (present(error)) then
+          error = error_l
+      else if (error_l .ne. ELPA_OK) then
+        write(error_unit,'(a)') "ELPA: Error in transform_generalized() and you did not check for errors!"
+      endif
 
       call self%get("solver", solver)
       if (solver .eq. ELPA_SOLVER_1STAGE) then
@@ -1257,6 +1266,13 @@ module elpa_impl
       else if (.not. success_l) then
         write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
       endif
+
+      call self%elpa_transform_back_generalized_d(b, q, sc_desc, error_l)
+      if (present(error)) then
+          error = error_l
+      else if (error_l .ne. ELPA_OK) then
+        write(error_unit,'(a)') "ELPA: Error in transform_back_generalized() and you did not check for errors!"
+      endif
     end subroutine
 
     !c> void elpa_generalized_eigenvectors_d(elpa_t handle, double *a, double *ev, double *q, int *error);
@@ -1321,10 +1337,18 @@ module elpa_impl
       integer             :: sc_desc(SC_DESC_LEN)
 
       integer, optional   :: error
+      integer             :: error_l
       integer(kind=c_int) :: solver
 #ifdef WANT_SINGLE_PRECISION_REAL
       logical             :: success_l
 
+      call self%elpa_transform_generalized_f(a, b, sc_desc, error_l)
+      if (present(error)) then
+          error = error_l
+      else if (error_l .ne. ELPA_OK) then
+        write(error_unit,'(a)') "ELPA: Error in transform_generalized() and you did not check for errors!"
+      endif
+
       call self%get("solver",solver)
       if (solver .eq. ELPA_SOLVER_1STAGE) then
         success_l = elpa_solve_evp_real_1stage_single_impl(self, a, ev, q)
@@ -1345,6 +1369,13 @@ module elpa_impl
       else if (.not. success_l) then
         write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
       endif
+
+      call self%elpa_transform_back_generalized_f(b, q, sc_desc, error_l)
+      if (present(error)) then
+          error = error_l
+      else if (error_l .ne. ELPA_OK) then
+        write(error_unit,'(a)') "ELPA: Error in transform_back_generalized() and you did not check for errors!"
+      endif
 #else
       print *,"This installation of the ELPA library has not been build with single-precision support"
       error = ELPA_ERROR
@@ -1415,9 +1446,17 @@ module elpa_impl
       integer                        :: sc_desc(SC_DESC_LEN)
 
       integer, optional              :: error
+      integer                        :: error_l
       integer(kind=c_int)            :: solver
       logical                        :: success_l
 
+      call self%elpa_transform_generalized_dc(a, b, sc_desc, error_l)
+      if (present(error)) then
+          error = error_l
+      else if (error_l .ne. ELPA_OK) then
+        write(error_unit,'(a)') "ELPA: Error in transform_generalized() and you did not check for errors!"
+      endif
+
       call self%get("solver", solver)
       if (solver .eq. ELPA_SOLVER_1STAGE) then
         success_l = elpa_solve_evp_complex_1stage_double_impl(self, a, ev, q)
@@ -1438,6 +1477,13 @@ module elpa_impl
       else if (.not. success_l) then
         write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
       endif
+
+      call self%elpa_transform_back_generalized_dc(b, q, sc_desc, error_l)
+      if (present(error)) then
+          error = error_l
+      else if (error_l .ne. ELPA_OK) then
+        write(error_unit,'(a)') "ELPA: Error in transform_back_generalized() and you did not check for errors!"
+      endif
     end subroutine
 
 
@@ -1505,10 +1551,18 @@ module elpa_impl
       integer                       :: sc_desc(SC_DESC_LEN)
 
       integer, optional             :: error
+      integer                       :: error_l
       integer(kind=c_int)           :: solver
 #ifdef WANT_SINGLE_PRECISION_COMPLEX
       logical                       :: success_l
 
+      call self%elpa_transform_generalized_fc(a, b, sc_desc, error_l)
+      if (present(error)) then
+          error = error_l
+      else if (error_l .ne. ELPA_OK) then
+        write(error_unit,'(a)') "ELPA: Error in transform_generalized() and you did not check for errors!"
+      endif
+
       call self%get("solver", solver)
       if (solver .eq. ELPA_SOLVER_1STAGE) then
         success_l = elpa_solve_evp_complex_1stage_single_impl(self, a, ev, q)
@@ -1529,6 +1583,13 @@ module elpa_impl
       else if (.not. success_l) then
         write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!"
       endif
+
+      call self%elpa_transform_back_generalized_fc(b, q, sc_desc, error_l)
+      if (present(error)) then
+          error = error_l
+      else if (error_l .ne. ELPA_OK) then
+        write(error_unit,'(a)') "ELPA: Error in transform_back_generalized() and you did not check for errors!"
+      endif
 #else
       print *,"This installation of the ELPA library has not been build with single-precision support"
       error = ELPA_ERROR

src/elpa_impl_template.F90

@@ -55,3 +55,36 @@
 
     end subroutine
 
+
+    subroutine elpa_transform_back_generalized_&
+            &ELPA_IMPL_SUFFIX&
+            &(self, b, q, sc_desc, error)
+        implicit none
+#include "general/precision_kinds.F90"
+        class(elpa_impl_t)  :: self
+#ifdef USE_ASSUMED_SIZE
+      MATH_DATATYPE(kind=rck) :: b(self%local_nrows, *), q(self%local_nrows, *)
+#else
+      MATH_DATATYPE(kind=rck) :: b(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols)
+#endif
+     integer                :: error
+     integer                :: sc_desc(9)
+
+     ! local helper array. TODO: do we want it this way? (do we need it? )
+     MATH_DATATYPE(kind=rck) :: tmp(self%local_nrows, self%local_ncols)
+
+     !todo: part of eigenvectors only
+#ifdef WITH_MPI
+     ! Q <= inv(U) * Q
+     call p&
+         &BLAS_CHAR&
+         &trmm("L", "U", "N", "N", self%na, self%na, &
+               ONE, b, 1, 1, sc_desc,  q, 1, 1, sc_desc)
+#else
+     call BLAS_CHAR&
+         &trmm("L", "U", "N", "N", self%na, self%na, &
+               ONE, b, self%na, q, self%na)
+#endif
+
+    end subroutine
+

test/Fortran/test.F90

@@ -344,8 +344,8 @@ program test
      !call prepare_matrix_random(na, myid, sc_desc, b, z, bs)
      ! TODO create random SPD matrix
      !diagonalElement = (2.546_rk, 0.0_rk)
-     diagonalElement = ONE
-     subdiagonalElement =  (0.0_rk, 0.0_rk)
+     diagonalElement = 2.546_rk * ONE
+     subdiagonalElement = ZERO
      call prepare_matrix_toeplitz(na, diagonalElement, subdiagonalElement, &
                                   d, sd, ds, sds, b, bs, nblk, np_rows, &
                                   np_cols, my_prow, my_pcol)
@@ -518,7 +518,12 @@ program test
 !     status = check_correctness_evp_numeric_residuals(na, nev, as, z, ev, sc_desc, nblk, myid, np_rows,np_cols, my_prow, my_pcol)
 !#elif defined(TEST_MATRIX_RANDOM)
      if (nev .ge. 1) then
+#if defined(TEST_GENERALIZED_EIGENPROBLEM)
+       status = check_correctness_evp_numeric_residuals(na, nev, as, z, ev, sc_desc, nblk, myid, np_rows,np_cols, my_prow, &
+                                                        my_pcol, bs)
+#else
        status = check_correctness_evp_numeric_residuals(na, nev, as, z, ev, sc_desc, nblk, myid, np_rows,np_cols, my_prow, my_pcol)
+#endif
      else
        ! zero eigenvectors and no analytic test => toeplitz
        status = check_correctness_eigenvalues_toeplitz(na, diagonalElement, &

test/shared/test_check_correctness_template.F90

@@ -45,12 +45,13 @@
     &MATH_DATATYPE&
     &_&
     &PRECISION&
-    & (na, nev, as, z, ev, sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol) result(status)
+    & (na, nev, as, z, ev, sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol, bs) result(status)
       implicit none
 #include "../../src/general/precision_kinds.F90"
       integer(kind=ik)                 :: status
       integer(kind=ik), intent(in)     :: na, nev, nblk, myid, np_rows, np_cols, my_prow, my_pcol
-      MATH_DATATYPE(kind=rck), intent(in)     :: as(:,:), z(:,:)
+      MATH_DATATYPE(kind=rck), intent(in)           :: as(:,:), z(:,:)
+      MATH_DATATYPE(kind=rck), intent(in), optional :: bs(:,:)
       real(kind=rk)                 :: ev(:)
       MATH_DATATYPE(kind=rck), dimension(size(as,dim=1),size(as,dim=2)) :: tmp1, tmp2
       MATH_DATATYPE(kind=rck)                :: xc
@@ -76,10 +77,13 @@
       real(kind=rk), parameter       :: tol_res_real_single      = 3e-3_rk
       real(kind=rk), parameter       :: tol_res_complex_double   = 5e-12_rk
       real(kind=rk), parameter       :: tol_res_complex_single   = 3e-3_rk
-      real(kind=rk), parameter       :: tol_res                  = tol_res_&
+      real(kind=rk)                  :: tol_res                  = tol_res_&
                                                                           &MATH_DATATYPE&
                                                                           &_&
                                                                           &PRECISION
+      ! precision of generalized problem is lower
+      real(kind=rk), parameter       :: generalized_penalty = 10.0_rk
+
       ! tolerance for the orthogonality test for different math type/precision setups
       real(kind=rk), parameter       :: tol_orth_real_double     = 5e-12_rk
       real(kind=rk), parameter       :: tol_orth_real_single     = 9e-4_rk
@@ -90,34 +94,50 @@
                                                                           &_&
                                                                           &PRECISION
 
+      if(present(bs)) then
+          tol_res = generalized_penalty * tol_res
+      endif
       status = 0
 
       ! 1. Residual (maximum of || A*Zi - Zi*EVi ||)
-      ! tmp1 =  A * Z
-      ! as is original stored matrix, Z are the EVs
 
-#ifdef WITH_MPI
-      call scal_PRECISION_GEMM('N', 'N', na, nev, na, ONE, as, 1, 1, sc_desc, &
-                  z, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
-#else /* WITH_MPI */
-      call PRECISION_GEMM('N','N',na,nev,na,ONE,as,na,z,na,ZERO,tmp1,na)
-#endif /* WITH_MPI */
-
-
-      ! tmp2 = Zi*EVi
-      tmp2(:,:) = z(:,:)
+      ! tmp1 = Zi*EVi
+      tmp1(:,:) = z(:,:)
       do i=1,nev
         xc = ev(i)
 #ifdef WITH_MPI
         call p&
             &BLAS_CHAR&
-            &scal(na, xc, tmp2, 1, i, sc_desc, 1)
+            &scal(na, xc, tmp1, 1, i, sc_desc, 1)
 #else /* WITH_MPI */
         call BLAS_CHAR&
-            &scal(na,xc,tmp2(:,i),1)
+            &scal(na,xc,tmp1(:,i),1)
 #endif /* WITH_MPI */
       enddo
 
+      ! for generalized EV problem, multiply by bs as well
+      ! tmp2 = B * tmp1
+      if(present(bs)) then
+#ifdef WITH_MPI
+      call scal_PRECISION_GEMM('N', 'N', na, nev, na, ONE, bs, 1, 1, sc_desc, &
+                  tmp1, 1, 1, sc_desc, ZERO, tmp2, 1, 1, sc_desc)
+#else /* WITH_MPI */
+      call PRECISION_GEMM('N','N',na,nev,na,ONE,bs,na,tmp1,na,ZERO,tmp2,na)
+#endif /* WITH_MPI */
+      else
+        ! normal eigenvalue problem .. no need to multiply
+        tmp2(:,:) = tmp1(:,:)
+      end if
+
+      ! tmp1 =  A * Z
+      ! as is original stored matrix, Z are the EVs
+#ifdef WITH_MPI
+      call scal_PRECISION_GEMM('N', 'N', na, nev, na, ONE, as, 1, 1, sc_desc, &
+                  z, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
+#else /* WITH_MPI */
+      call PRECISION_GEMM('N','N',na,nev,na,ONE,as,na,z,na,ZERO,tmp1,na)
+#endif /* WITH_MPI */
+
       !  tmp1 = A*Zi - Zi*EVi
       tmp1(:,:) =  tmp1(:,:) - tmp2(:,:)
 
@@ -166,7 +186,9 @@
       endif
 
       ! 2. Eigenvector orthogonality
-
+      !TODO for the generalized eigenvector problem, the orthogonality test has to be altered
+      !TODO at the moment, it is skipped
+      if(.not. present(bs)) then
       ! tmp1 = Z**T * Z
       tmp1 = 0
 #ifdef WITH_MPI
@@ -223,6 +245,8 @@
           status = 1
         endif
       endif
+
+      endif  ! skiping test of orthogonality for generalized eigenproblem
     end function