text: /******************************************************************************
  * Copyright (C) 2021, Advanced Micro Devices, Inc. All rights reserved.
  *******************************************************************************/
  
  /*! @file libflame_hbgv.cc
   *  @brief Test application to validate hbgv() using CPP template interface.
   *  */
  
  #include <gtest/gtest.h>
  #include "main.h"
  #include "libflame_test.hh"
  
  /*! @brief  hbgv_test is function template for hbgv() functions.
              T can be scomplex, dcomplex
              Ta can be float, double.
   * @details
   * \b Purpose:
      \verbatim
  	  hbgv_test is function template for hbgv() functions.
  	  T can be scomplex, dcomplex
      Ta can be float, double.
  	  
  	  hbgv_test() function template calls C and CPP based lbrary APIs with
  	  valid test values, calculate the differences in output if INFO is >= 0.
      And passses the test case if difference is <= threshold.
      Fails the test case if difference > threshold or INFO < 0.
      
      Complex reference:
  	  http://www.netlib.org/lapack/explore-html/d9/d98/group__complex_o_t_h_e_reigen_gae30c26efa0a7b94048c00cad17532044.html#gae30c26efa0a7b94048c00cad17532044
  	  Complex double reference:
  	  http://www.netlib.org/lapack/explore-html/db/d61/group__complex16_o_t_h_e_reigen_ga76a8cfc758f8dc17ac37f6eed2ef18a4.html#ga76a8cfc758f8dc17ac37f6eed2ef18a4
      \endverbatim
  	
   * @param[in] IP
            IP is INTEGER
            Used to pass Index of Eigen Parameters array present in config file.
  
   * @return VOID
             Nothing.
   * */
  template< typename T, typename Ta >
  void hbgv_test(int ip)
  {
    typedef integer (*fptr_NL_LAPACK_hbgv)(char* jobz, char* uplo, integer* n,
                        integer* ka, integer* kb, T* ab, integer* ldab, T* bb,
                        integer* ldbb, Ta* w, T* z, integer* ldz, T* work,
                        Ta* rwork, integer* info);
    fptr_NL_LAPACK_hbgv hbgv_ref = NULL;
    
    // Initialise random number generators with timestamp.
    srand (SRAND_SEED_VALUE);
    
    /* N is INTEGER
            The order of the matrices A and B.  N >= 0.*/
    integer n = eig_paramslist[ip].n;
    if (n < 0) {
      PRINTF("n < 0 but should be: n >= 0. Please correct the input data.\n");
    }
    
    /* JOBZ is CHARACTER*1
            = 'N':  Compute eigenvalues only;
            = 'V':  Compute eigenvalues and eigenvectors.*/
    char jobz = eig_paramslist[ip].jobz;
    if ((jobz != 'N') && (jobz != 'V')) {
      PRINTF("jobz should be N or V. Please correct the input data.\n");
    }
    
    /* UPLO is CHARACTER*1
            = 'U':  Upper triangles of A and B are stored;
            = 'L':  Lower triangles of A and B are stored.*/
    char uplo = eig_paramslist[ip].uplo;
    if ((uplo != 'U') && (uplo != 'L')) {
      PRINTF("jobz should be N or V. Please correct the input data.\n");
    }
  
    /* KA is INTEGER
            The number of superdiagonals of the matrix A if UPLO = 'U',
            or the number of subdiagonals if UPLO = 'L'. KA >= 0.*/
    integer ka = 0;
    if (uplo == 'U') {
  	  ka = eig_paramslist[ip].sda;
    } else if (uplo == 'L') {
  	  ka = eig_paramslist[ip].subda;
    }
    if (ka < 0) {
      PRINTF("ka < 0 but should be: KA >= 0. Please correct the input data.\n");
    }
    
    /* KB is INTEGER
            The number of superdiagonals of the matrix B if UPLO = 'U',
            or the number of subdiagonals if UPLO = 'L'. KB >= 0.*/
    integer kb = 0;
    if (uplo == 'U') {
  	  kb = eig_paramslist[ip].sdb;
    } else if (uplo == 'L') {
  	  kb = eig_paramslist[ip].subdb;
    }
    if (kb < 0) {
      PRINTF("kb < 0 but should be: KB >= 0. Please correct the input data.\n");
    }
    
    /* LDAB is INTEGER
            The leading dimension of the array AB.  LDAB >= KA+1.*/
    integer ldab = eig_paramslist[ip].ldab_hbgv;
    if (ldab < (ka + 1)) {
      PRINTF("ldab < (ka + 1) but should be: LDAB >= (KA + 1). Please correct" \
             " the input data.\n");
    }
    
    // AB is COMPLEX or COMPLEX*16 array, dimension (LDAB,N)
    T *abbuff, *abrefbuff;
    allocate_init_buffer(abbuff, abrefbuff, ldab * n);
    
    /* LDBB is INTEGER
            The leading dimension of the array BB.  LDBB >= KB+1.*/
    integer ldbb = eig_paramslist[ip].ldbb;
    if (ldbb < (kb + 1)) {
      PRINTF("ldbb < (kb + 1) but should be: LDBB >= (KB + 1). Please correct" \
             " the input data.\n");
    }
    
    /* BB is COMPLEX or COMPLEX*16 array, dimension (LDBB,N) */
    T *bbbuff, *bbrefbuff;
    allocate_init_buffer(bbbuff, bbrefbuff, ldbb * n);
    
    Ta *wbuff, *wrefbuff;
    // W is REAL or DOUBLE PRECISION array, dimension (N)
    allocate_init_buffer(wbuff, wrefbuff, n, 0);
    
    /* LDZ is INTEGER
            The leading dimension of the array Z.  LDZ >= 1, and if
            JOBZ = 'V', LDZ >= N.*/
    integer ldz = eig_paramslist[ip].ldz_hbgv;
    if (ldz < 1) {
      PRINTF("ldz < 1 but should be: LDZ >= 1. Please correct the input data.\n");
    }
    if ((jobz == 'V') && (ldz < n)) {
      PRINTF("If jobz = V and ldz < n but should be: LDZ >= N. Please correct" \
             " the input data.\n");
    }
    
    /* Z is COMPLEX or COMPLEX*16 array, dimension (LDZ, N)
            If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
            eigenvectors, with the i-th column of Z holding the
            eigenvector associated with W(i). The eigenvectors are
            normalized so that Z**H*B*Z = I.
            If JOBZ = 'N', then Z is not referenced.*/
    T *zbuff, *zrefbuff; // output buffer
    allocate_init_buffer(zbuff, zrefbuff, ldz * n, 0);
    
    // WORK is COMPLEX or COMPLEX*16  array, dimension (N)
    T *workbuff = NULL, *workrefbuff = NULL;
    allocate_init_buffer(workbuff, workrefbuff, n, 0);
    
    // RWORK is REAL or DOUBLE PRECISION array, dimension (LRWORK)
    Ta *rworkbuff = NULL, *rworkrefbuff = NULL;
    allocate_init_buffer(rworkbuff, rworkrefbuff, 3 * n, 0);
    
    // Print input values other than arrays.
    #if (defined(PRINT_INPUT_VALUES) && (PRINT_INPUT_VALUES == 1))
      PRINTF("\nPrinting all Input values other than array contents...\n");
      PRINTF("jobz = %c\n", jobz);
      PRINTF("uplo = %c\n", uplo);
      PRINTF("n = %d\n", n);
      PRINTF("ka = %d\n", ka);
      PRINTF("kb = %d\n", kb);
      PRINTF("ldab = %d\n", ldab);
      PRINTF("Size of AB array (ldab*n) = %d\n", ldab*n);
      PRINTF("ldbb = %d\n", ldbb);
      PRINTF("Size of BB array (ldbb*n) = %d\n", ldbb*n);
      PRINTF("ldz = %d\n", ldz);
      PRINTF("Size of Z array (ldz*n) = %d\n", ldz*n);
      PRINTF("Size of WORK array (n) = %d\n", n);
      PRINTF("Size of RWORK array (3*n) = %d\n", 3*n);
    #endif
  
    #if (defined(PRINT_ARRAYS) && (PRINT_ARRAYS == 1))
    // Array to store array name to print.
    char arrayname[20] = "";
    integer arraysize = sizeof(arrayname);
    #endif
    
    #if (defined(PRINT_ARRAYS) && (PRINT_ARRAYS == 1) && \
        defined(PRINT_INPUT_ARRAYS) && (PRINT_INPUT_ARRAYS == 1))
      // Print all input arrays if PRINT_INPUT_ARRAYS macro is enabled
      PRINTF("\nPrinting all Input arrays contents...\n");
      
      // Prints AB array contents
      strncpy(arrayname, "AB input", arraysize);
      print_array<T>(arrayname, abbuff, ldab * n);
      strncpy(arrayname, "AB ref input", arraysize);
      print_array<T>(arrayname, abrefbuff, ldab * n);
      
      // Prints BB array contents
      strncpy(arrayname, "BB input", arraysize);
      print_array<T>(arrayname, bbbuff, ldbb * n);
      strncpy(arrayname, "BB ref input", arraysize);
      print_array<T>(arrayname, bbrefbuff, ldbb * n);
      
      // Prints W array contents
      strncpy(arrayname, "W input", arraysize);
      print_array<Ta>(arrayname, wbuff, n);
      strncpy(arrayname, "W ref input", arraysize);
      print_array<Ta>(arrayname, wrefbuff, n);
      
      // Prints Z array contents
      strncpy(arrayname, "Z input", arraysize);
      print_array<T>(arrayname, zbuff, ldz * n);
      strncpy(arrayname, "Z ref input", arraysize);
      print_array<T>(arrayname, zrefbuff, ldz * n);
      
      // Prints WORK array contents
      strncpy(arrayname, "WORK input", arraysize);
      print_array<T>(arrayname, workbuff, n);
      strncpy(arrayname, "WORK ref input", arraysize);
      print_array<T>(arrayname, workrefbuff, n);
      
      // Prints RWORK array contents
      strncpy(arrayname, "RWORK input", arraysize);
      print_array<Ta>(arrayname, rworkbuff, 3 * n);
      strncpy(arrayname, "RWORK ref input", arraysize);
      print_array<Ta>(arrayname, rworkrefbuff, 3 * n);
    #endif
    
    integer info_cpp = -1;
    // Call CPP function
    libflame::hbgv<T, Ta>(&jobz, &uplo, &n, &ka, &kb, abbuff, &ldab, bbbuff,
                &ldbb, wbuff, zbuff, &ldz, workbuff, rworkbuff, &info_cpp);
  
    // Call C function
    /* Check the typename T passed to this function template and call respective
       function.*/
    if (typeid(T) == typeid(scomplex)) {
      hbgv_ref = (fptr_NL_LAPACK_hbgv)dlsym(lapackModule, "chbgv_");
    } else if (typeid(T) == typeid(dcomplex)) {
      hbgv_ref = (fptr_NL_LAPACK_hbgv)dlsym(lapackModule, "zhbgv_");
    } else {
  	  PRINTF("Invalid typename is passed to %s() function template.\n",
             __FUNCTION__);
    }
    
    if (hbgv_ref == NULL) {
      PRINTF("Could not get the symbol. Exiting...\n");
  	  closelibs();
      exit(-1);
    }
    integer info_ref = -1;
    
    hbgv_ref(&jobz, &uplo, &n, &ka, &kb, abrefbuff, &ldab, bbrefbuff, &ldbb,
          wrefbuff, zrefbuff, &ldz, workrefbuff, rworkrefbuff, &info_ref);
    PRINTF ("info_cpp: %d, info_ref: %d\n", info_cpp, info_ref);
    
    // Calculate the differences of buffers.
    if ((info_cpp >= 0) && (info_ref >= 0)) {
      #if (defined(PRINT_ARRAYS) && (PRINT_ARRAYS == 1) && \
          defined(PRINT_OUTPUT_ARRAYS) && (PRINT_OUTPUT_ARRAYS == 1))
        // Print all output arrays if PRINT_OUTPUT_ARRAYS macro is enabled
        PRINTF("\nPrinting all Output arrays contents...\n");
        
        // Prints AB array contents
        strncpy(arrayname, "AB output", arraysize);
        print_array<T>(arrayname, abbuff, ldab * n);
        strncpy(arrayname, "AB ref output", arraysize);
        print_array<T>(arrayname, abrefbuff, ldab * n);
        
        // Prints BB array contents
        strncpy(arrayname, "BB output", arraysize);
        print_array<T>(arrayname, bbbuff, ldbb * n);
        strncpy(arrayname, "BB ref output", arraysize);
        print_array<T>(arrayname, bbrefbuff, ldbb * n);
        
        // Prints W array contents
        strncpy(arrayname, "W output", arraysize);
        print_array<Ta>(arrayname, wbuff, n);
        strncpy(arrayname, "W ref output", arraysize);
        print_array<Ta>(arrayname, wrefbuff, n);
        
        // Prints Z array contents
        strncpy(arrayname, "Z output", arraysize);
        print_array<T>(arrayname, zbuff, ldz * n);
        strncpy(arrayname, "Z ref output", arraysize);
        print_array<T>(arrayname, zrefbuff, ldz * n);
        
        // Prints WORK array contents
        strncpy(arrayname, "WORK output", arraysize);
        print_array<T>(arrayname, workbuff, n);
        strncpy(arrayname, "WORK ref output", arraysize);
        print_array<T>(arrayname, workrefbuff, n);
        
        // Prints RWORK array contents
        strncpy(arrayname, "RWORK output", arraysize);
        print_array<Ta>(arrayname, rworkbuff, 3 * n);
        strncpy(arrayname, "RWORK ref output", arraysize);
        print_array<Ta>(arrayname, rworkrefbuff, 3 * n);
      #endif
      
      double diff = computeError<T>(ldab, n, abrefbuff, abbuff);
      diff += computeError<T>(ldbb, n, bbrefbuff, bbbuff);
      if (jobz == 'V') {
        diff += computeError<T>(ldz, n, zbuff, zrefbuff);
      }
      diff += computeError<Ta>(1, n, wbuff, wrefbuff);
      diff += computeError<T>(1, n, workbuff, workrefbuff);
      diff += computeError<Ta>(3, n, rworkbuff, rworkrefbuff);
      PRINTF("diff: %lf\n", diff);
      EXPECT_NEAR(0.0, abs(diff), SYM_EIGEN_THRESHOLD);
    } else {
      PRINTF("Info returned by CPP or C API is not successful to compare" \
              " differences.\n");
      EXPECT_FALSE((info_cpp < 0) || (info_ref < 0));
    }
    
    // Free up the buffers.
    delete[] abbuff; delete[] abrefbuff;
    delete[] bbbuff; delete[] bbrefbuff;
    delete[] wbuff; delete[] wrefbuff;
    delete[] zbuff; delete[] zrefbuff;
    delete[] workbuff; delete[] workrefbuff;
    delete[] rworkbuff; delete[] rworkrefbuff;
  }
  
  /* Use TEST macro and call C++ test function template with
     scomplex and float as typenames.*/
  TEST(LAPACKCPP_hbgv, CHBGV) {
    for (short int index = 0; index < NUM_SUB_TESTS; index++) {
  	  PRINTF("index: %d\n", index);
      hbgv_test<scomplex, float> (index);
    }
  }
  
  /* Use TEST macro and call C++ test function template with
     dcomplex and double as typenames.*/
  TEST(LAPACKCPP_hbgv, ZHBGV) {
    for (short int index = 0; index < NUM_SUB_TESTS; index++) {
  	  PRINTF("index: %d\n", index);
      hbgv_test<dcomplex, double> (index);
    }
  }