!    This file is part of ELPA.
!
!    The ELPA library was originally created by the ELPA consortium,
!    consisting of the following organizations:
!
!    - Max Planck Computing and Data Facility (MPCDF), formerly known as
!      Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
!    - Bergische Universität Wuppertal, Lehrstuhl für angewandte
!      Informatik,
!    - Technische Universität München, Lehrstuhl für Informatik mit
!      Schwerpunkt Wissenschaftliches Rechnen ,
!    - Fritz-Haber-Institut, Berlin, Abt. Theorie,
!    - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
!      Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
!      and
!    - IBM Deutschland GmbH
!
!
!    More information can be found here:
!    http://elpa.mpcdf.mpg.de/
!
!    ELPA is free software: you can redistribute it and/or modify
!    it under the terms of the version 3 of the license of the
!    GNU Lesser General Public License as published by the Free
!    Software Foundation.
!
!    ELPA is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public License
!    along with ELPA.  If not, see <http://www.gnu.org/licenses/>
!
!    ELPA reflects a substantial effort on the part of the original
!    ELPA consortium, and we ask you to respect the spirit of the
!    license that we chose: i.e., please contribute any changes you
!    may have back to the original ELPA library distribution, and keep
!    any derivatives of ELPA under the same license that we chose for
!    the original distribution, the GNU Lesser General Public License.
!
! Author: A. Marek, MPCDF

    function check_correctness_&
    &MATH_DATATYPE&
    &_&
    &PRECISION&
    & (na, nev, as, z, ev, sc_desc, myid) result(status)
      implicit none
      integer(kind=ik)                 :: status
      integer(kind=ik), intent(in)     :: na, nev, myid
#if REALCASE == 1
      real(kind=C_DATATYPE_KIND), intent(in)     :: as(:,:), z(:,:)
      real(kind=C_DATATYPE_KIND)                 :: ev(:)
      real(kind=C_DATATYPE_KIND), dimension(size(as,dim=1),size(as,dim=2)) :: tmp1, tmp2

#ifdef DOUBLE_PRECISION_REAL
      real(kind=C_DATATYPE_KIND), parameter      :: ZERO=0.0_rk8, ONE = 1.0_rk8
#else
      real(kind=C_DATATYPE_KIND), parameter      :: ZERO=0.0_rk4, ONE = 1.0_rk4
#endif

#ifndef WITH_MPI

#ifdef DOUBLE_PRECISION_REAL
      real(kind=C_DATATYPE_KIND)                 :: dnrm2
#else
      real(kind=C_DATATYPE_KIND)                 :: snrm2
#endif

#endif
#endif /* REALCASE */

#if COMPLEXCASE == 1
      complex(kind=C_DATATYPE_KIND), intent(in)    :: as(:,:), z(:,:)
      real(kind=C_DATATYPE_KIND)                   :: ev(:)
      complex(kind=C_DATATYPE_KIND), dimension(size(as,dim=1),size(as,dim=2)) :: tmp1, tmp2
      complex(kind=C_DATATYPE_KIND)                :: xc
#ifdef DOUBLE_PRECISION_COMPLEX
      complex(kind=C_DATATYPE_KIND), parameter     :: ZERO = (0.0_rk8,0.0_rk8), ONE = (1.0_rk8,0.0_rk8)
#ifndef WITH_MPI
      complex(kind=C_DATATYPE_KIND)                :: zdotc, cdotc
#endif

#else /* DOUBLE_PRECISION_COMPLEX */
      complex(kind=C_DATATYPE_KIND), parameter     :: ZERO = (0.0_rk4,0.0_rk4), ONE = (1.0_rk4,0.0_rk4)
#ifndef WITH_MPI
      complex(kind=C_DATATYPE_KIND)                :: zdotc, cdotc
#endif

#endif /* DOUBLE_PRECISION_COMPLEX */

#endif /* COMPLEXCASE */

      integer(kind=ik)                 :: sc_desc(:)

      integer(kind=ik)                 :: i
      real(kind=C_DATATYPE_KIND)                   :: err, errmax

      integer :: mpierr

      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(:,:)
      do i=1,nev
#if REALCASE == 1
#ifdef WITH_MPI

#ifdef DOUBLE_PRECISION_REAL
        call pdscal(na, ev(i), tmp2, 1, i, sc_desc, 1)
#else
        call psscal(na, ev(i), tmp2, 1, i, sc_desc, 1)
#endif

#else /* WITH_MPI */

#ifdef DOUBLE_PRECISION_REAL
        call dscal(na,ev(i),tmp2(:,i),1)
#else
        call sscal(na,ev(i),tmp2(:,i),1)
#endif

#endif /* WITH_MPI */
#endif /* REALCASE */

#if COMPLEXCASE == 1
        xc = ev(i)
#ifdef WITH_MPI

#ifdef DOUBLE_PRECISION_COMPLEX
        call pzscal(na, xc, tmp2, 1, i, sc_desc, 1)
#else
        call pcscal(na, xc, tmp2, 1, i, sc_desc, 1)
#endif

#else /* WITH_MPI */

#ifdef DOUBLE_PRECISION_COMPLEX
        call zscal(na,xc,tmp2(1,i),1)
#else
        call cscal(na,xc,tmp2(1,i),1)
#endif

#endif /* WITH_MPI */
#endif /* COMPLEXCASE */
      enddo

      !  tmp1 = A*Zi - Zi*EVi
      tmp1(:,:) =  tmp1(:,:) - tmp2(:,:)

      ! Get maximum norm of columns of tmp1
#if REALCASE == 1
#ifdef DOUBLE_PRECISION_REAL
      errmax = 0.0_rk8
#else
      errmax = 0.0_rk4
#endif
#endif
#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
      errmax = 0.0_rk8
#else
      errmax = 0.0_rk4
#endif
#endif

      do i=1,nev
#if REALCASE == 1
#ifdef DOUBLE_PRECISION_REAL
        err = 0.0_rk8
#else
        err = 0.0_rk4
#endif

#ifdef WITH_MPI

#ifdef DOUBLE_PRECISION_REAL
        call pdnrm2(na, err, tmp1, 1, i, sc_desc, 1)
#else
        call psnrm2(na, err, tmp1, 1, i, sc_desc, 1)
#endif

#else /* WITH_MPI */

#ifdef DOUBLE_PRECISION_REAL
!        call dnrm2(na,err,tmp1,1,i,sc_desc,1)
        err = dnrm2(na,tmp1(1,i),1)
#else
        err = snrm2(na,tmp1(1,i),1)
#endif

#endif /* WITH_MPI */
        errmax = max(errmax, err)
#endif /* REALCASE */

#if COMPLEXCASE == 1
        xc = 0
#ifdef WITH_MPI

#ifdef DOUBLE_PRECISION_COMPLEX
        call pzdotc(na, xc, tmp1, 1, i, sc_desc, 1, tmp1, 1, i, sc_desc, 1)
#else
        call pcdotc(na, xc, tmp1, 1, i, sc_desc, 1, tmp1, 1, i, sc_desc, 1)
#endif

#else /* WITH_MPI */

#ifdef DOUBLE_PRECISION_COMPLEX
        xc = zdotc(na,tmp1,1,tmp1,1)
#else
        xc = cdotc(na,tmp1,1,tmp1,1)
#endif

#endif /* WITH_MPI */


#ifdef DOUBLE_PRECISION_COMPLEX
        errmax = max(errmax, sqrt(real(xc,kind=rk8)))
#else
        errmax = max(errmax, sqrt(real(xc,kind=rk4)))
#endif
#endif /* COMPLEXCASE */
      enddo

      ! Get maximum error norm over all processors
      err = errmax
#ifdef WITH_MPI
      call mpi_allreduce(err, errmax, 1, MPI_REAL_PRECISION, MPI_MAX, MPI_COMM_WORLD, mpierr)
#else /* WITH_MPI */
      errmax = err
#endif /* WITH_MPI */

      if (myid==0) print *,'Error Residual     :',errmax
#if REALCASE == 1
#ifdef DOUBLE_PRECISION_REAL
      if (errmax .gt. 5e-12_rk8 .or. errmax .eq. 0.0_rk8) then
#else
      if (errmax .gt. 3e-3_rk4 .or. errmax .eq. 0.0_rk4) then
#endif
#endif
#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
      if (errmax .gt. 5e-12_rk8 .or. errmax .eq. 0.0_rk8) then
#else
      if (errmax .gt. 3e-3_rk4 .or. errmax .eq. 0.0_rk4) then
#endif
#endif

        status = 1
      endif

      ! 2. Eigenvector orthogonality

      ! tmp1 = Z**T * Z
      tmp1 = 0
#ifdef WITH_MPI
      call scal_PRECISION_GEMM(BLAS_TRANS_OR_CONJ, 'N', nev, nev, na, ONE, z, 1, 1, &
                        sc_desc, z, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
#else /* WITH_MPI */
      call PRECISION_GEMM(BLAS_TRANS_OR_CONJ,'N',nev,nev,na,ONE,z,na,z,na,ZERO,tmp1,na)
#endif /* WITH_MPI */

      ! Initialize tmp2 to unit matrix
      tmp2 = 0
#if REALCASE == 1

#ifdef WITH_MPI

#ifdef DOUBLE_PRECISION_REAL
      call pdlaset('A', nev, nev, 0.0_rk8, 1.0_rk8, tmp2, 1, 1, sc_desc)
#else
      call pslaset('A', nev, nev, 0.0_rk4, 1.0_rk4, tmp2, 1, 1, sc_desc)
#endif

#else /* WITH_MPI */

#ifdef DOUBLE_PRECISION_REAL
      call dlaset('A',nev,nev,0.0_rk8,1.0_rk8,tmp2,na)
#else
      call slaset('A',nev,nev,0.0_rk4,1.0_rk4,tmp2,na)
#endif

#endif /* WITH_MPI */

#endif /* REALCASE */

#if COMPLEXCASE == 1
#ifdef WITH_MPI

#ifdef DOUBLE_PRECISION_COMPLEX
      call pzlaset('A', nev, nev, ZERO, ONE, tmp2, 1, 1, sc_desc)
#else
      call pclaset('A', nev, nev, ZERO, ONE, tmp2, 1, 1, sc_desc)
#endif

#else /* WITH_MPI */

#ifdef DOUBLE_PRECISION_COMPLEX
      call zlaset('A',nev,nev,ZERO,ONE,tmp2,na)
#else
      call claset('A',nev,nev,ZERO,ONE,tmp2,na)
#endif

#endif /* WITH_MPI */
#endif /* COMPLEXCASE */


      !      ! tmp1 = Z**T * Z - Unit Matrix
      tmp1(:,:) =  tmp1(:,:) - tmp2(:,:)

      ! Get maximum error (max abs value in tmp1)
      err = maxval(abs(tmp1))
#ifdef WITH_MPI
      call mpi_allreduce(err, errmax, 1, MPI_REAL_PRECISION, MPI_MAX, MPI_COMM_WORLD, mpierr)
#else /* WITH_MPI */
      errmax = err
#endif /* WITH_MPI */

      if (myid==0) print *,'Error Orthogonality:',errmax
#if REALCASE == 1
#ifdef DOUBLE_PRECISION_REAL
      if (errmax .gt. 5e-12_rk8 .or. errmax .eq. 0.0_rk8) then
#else
      if (errmax .gt. 9e-4_rk4 .or. errmax .eq. 0.0_rk4) then
#endif
#endif
#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
      if (errmax .gt. 5e-12_rk8 .or. errmax .eq. 0.0_rk8) then
#else
      if (errmax .gt. 9e-4_rk4 .or. errmax .eq. 0.0_rk4) then
#endif
#endif

        status = 1
      endif
    end function


#if REALCASE == 1

#ifdef DOUBLE_PRECISION_REAL
    !c> int check_correctness_real_double_f(int na, int nev, int na_rows, int na_cols,
    !c>                                         double *as, double *z, double *ev,
    !c>                                         int sc_desc[9], int myid);
#else
    !c> int check_correctness_real_single_f(int na, int nev, int na_rows, int na_cols,
    !c>                                         float *as, float *z, float *ev,
    !c>                                         int sc_desc[9], int myid);
#endif

#endif /* REALCASE */

#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
    !c> int check_correctness_complex_double_f(int na, int nev, int na_rows, int na_cols,
    !c>                                         complex double *as, complex double *z, double *ev,
    !c>                                         int sc_desc[9], int myid);
#else
    !c> int check_correctness_complex_single_f(int na, int nev, int na_rows, int na_cols,
    !c>                                         complex float *as, complex float *z, float *ev,
    !c>                                         int sc_desc[9], int myid);
#endif
#endif /* COMPLEXCASE */

function check_correctness_&
&MATH_DATATYPE&
&_&
&PRECISION&
&_f (na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid) result(status) &
      bind(C,name="check_correctness_&
      &MATH_DATATYPE&
      &_&
      &PRECISION&
      &_f")

      use iso_c_binding

      implicit none

      integer(kind=c_int)            :: status
      integer(kind=c_int), value     :: na, nev, myid, na_rows, na_cols
#if REALCASE == 1
      real(kind=C_DATATYPE_KIND)     :: as(1:na_rows,1:na_cols), z(1:na_rows,1:na_cols)
#endif

#if COMPLEXCASE == 1
      complex(kind=C_DATATYPE_KIND) :: as(1:na_rows,1:na_cols), z(1:na_rows,1:na_cols)
#endif
      real(kind=C_DATATYPE_KIND)    :: ev(1:na)
      integer(kind=c_int)            :: sc_desc(1:9)

      status = check_correctness_&
      &MATH_DATATYPE&
      &_&
      &PRECISION&
      & (na, nev, as, z, ev, sc_desc, myid)

    end function

    function check_correctness_eigenvalues_toeplitz_&
    &MATH_DATATYPE&
    &_&
    &PRECISION&
    & (na, diagonalElement, subdiagonalElement, ev, z, myid) result(status)
      use iso_c_binding
      implicit none

      integer               :: status, ii, j, myid
      integer, intent(in)   :: na
      real(kind=C_DATATYPE_KIND) :: diagonalElement, subdiagonalElement
      real(kind=C_DATATYPE_KIND) :: ev_analytic(na), ev(na)
#if REALCASE == 1
      real(kind=C_DATATYPE_KIND) :: z(:,:)
#endif
#if COMPLEXCASE == 1
      complex(kind=C_DATATYPE_KIND) :: z(:,:)
#endif

#if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
      real(kind=C_DATATYPE_KIND), parameter   :: pi = 3.141592653589793238462643383279_c_double
#else
      real(kind=C_DATATYPE_KIND), parameter   :: pi = 3.1415926535897932_c_float
#endif
      real(kind=C_DATATYPE_KIND)              :: tmp, maxerr
      integer                                 :: loctmp
      status = 0

     ! analytic solution
     do ii=1, na
#if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
       ev_analytic(ii) = diagonalElement + 2.0_c_double * &
                         subdiagonalElement *cos( pi*real(ii,kind=c_double)/ &
			 real(na+1,kind=c_double) )
#else
       ev_analytic(ii) = diagonalElement + 2.0_c_float * &
                         subdiagonalElement *cos( pi*real(ii,kind=c_float)/ &
                         real(na+1,kind=c_float) )
#endif
     enddo

     ! sort analytic solution:

     ! this hack is neither elegant, nor optimized: for huge matrixes it might be expensive
     ! a proper sorting algorithmus might be implemented here

     tmp    = minval(ev_analytic)
     loctmp = minloc(ev_analytic, 1)

     ev_analytic(loctmp) = ev_analytic(1)
     ev_analytic(1) = tmp
     do ii=2, na
       tmp = ev_analytic(ii)
       do j= ii, na
         if (ev_analytic(j) .lt. tmp) then
           tmp    = ev_analytic(j)
           loctmp = j
         endif
       enddo
       ev_analytic(loctmp) = ev_analytic(ii)
       ev_analytic(ii) = tmp
     enddo

     ! compute a simple error max of eigenvalues
     maxerr = 0.0
     maxerr = maxval( (ev(:) - ev_analytic(:))/ev_analytic(:) , 1)

#if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
     if (maxerr .gt. 8.e-13_c_double) then
#else
     if (maxerr .gt. 8.e-4_c_float) then
#endif
       status = 1
       if (myid .eq. 0) then
         print *,"Eigenvalues differ from analytic solution: maxerr = ",maxerr
       endif
     endif
    end function

! vim: syntax=fortran