/*  Accelerate: zheev and zhegv with UPLO = 'U' return wrong eigenvectors, with INFO = 0,
 *  for a Hermitian matrix of order N >= 130 whose last two (or more) rows and columns are
 *  exactly zero. The eigenvalues are correct. UPLO = 'L', zheevd and dsyev are correct on
 *  the same matrices.
 *
 *  Build and run, legacy LAPACK interface:
 *    clang -O2 zheev_upper_bug.c -o zheev_upper_bug -framework Accelerate && ./zheev_upper_bug
 *  Build and run, new LAPACK interface:
 *    clang -O2 -DACCELERATE_NEW_LAPACK zheev_upper_bug.c -o zheev_upper_bug_new -framework Accelerate && ./zheev_upper_bug_new
 *
 *  The matrix: random Hermitian, entries uniform in [-0.5,0.5), bandwidth 11 (a dense
 *  matrix fails as well), then the last NZERO rows and columns set to zero. The test
 *  reports max |V^H V - 1| and max |A V - V diag(W)| for the returned eigenvectors V.
 */
#include 
#include 
#include 
#include 
#include 

#ifdef ACCELERATE_NEW_LAPACK
typedef __LAPACK_int          lint;
typedef __LAPACK_double_complex zc;
#define RE(z) (((double*)&(z))[0])
#define IM(z) (((double*)&(z))[1])
static const char *interface_name = "new LAPACK interface (ACCELERATE_NEW_LAPACK)";
#else
typedef __CLPK_integer        lint;
typedef __CLPK_doublecomplex  zc;
#define RE(z) ((z).r)
#define IM(z) ((z).i)
static const char *interface_name = "legacy CLAPACK interface";
#endif

static unsigned long seed = 12345;
static double rnd(void)
{
  seed = (16807 * seed) % 2147483647UL;
  return (double)seed / 2147483647.0 - 0.5;
}

/*  column major, A[i + j*n] */
static void build(zc *A, int n, int band, int nzero)
{
  memset(A, 0, sizeof(zc) * n * n);
  for (int j = 0; j < n; j++)
  {
    RE(A[j + j*n]) = rnd();
    for (int i = 0; i < j; i++)
      if (j - i <= band)
      {
        double re = rnd(), im = rnd();
        RE(A[i + j*n]) = re; IM(A[i + j*n]) =  im;
        RE(A[j + i*n]) = re; IM(A[j + i*n]) = -im;
      }
  }
  for (int z = n - nzero; z < n; z++)
    for (int k = 0; k < n; k++)
    {
      RE(A[z + k*n]) = IM(A[z + k*n]) = 0;
      RE(A[k + z*n]) = IM(A[k + z*n]) = 0;
    }
}

/*  max |V^H V - 1| and max |A V - V diag(w)| */
static void check(const zc *A, const zc *V, const double *w, int n, double *orth, double *res)
{
  *orth = 0; *res = 0;
  for (int a = 0; a < n; a++)
    for (int b = 0; b < n; b++)
    {
      double re = 0, im = 0;
      for (int k = 0; k < n; k++)
      {
        re += RE(V[k + a*n]) * RE(V[k + b*n]) + IM(V[k + a*n]) * IM(V[k + b*n]);
        im += RE(V[k + a*n]) * IM(V[k + b*n]) - IM(V[k + a*n]) * RE(V[k + b*n]);
      }
      if (a == b) re -= 1;
      *orth = fmax(*orth, hypot(re, im));
    }
  for (int a = 0; a < n; a++)
    for (int i = 0; i < n; i++)
    {
      double re = 0, im = 0;
      for (int k = 0; k < n; k++)
      {
        re += RE(A[i + k*n]) * RE(V[k + a*n]) - IM(A[i + k*n]) * IM(V[k + a*n]);
        im += RE(A[i + k*n]) * IM(V[k + a*n]) + IM(A[i + k*n]) * RE(V[k + a*n]);
      }
      re -= w[a] * RE(V[i + a*n]);
      im -= w[a] * IM(V[i + a*n]);
      *res = fmax(*res, hypot(re, im));
    }
}

/*  routine: 0 zheev, 1 zhegv with B = identity, 2 zheevd */
static void run(int routine, char uplo, int n, int band, int nzero)
{
  zc *A = malloc(sizeof(zc) * n * n), *V = malloc(sizeof(zc) * n * n), *B = malloc(sizeof(zc) * n * n);
  double *w = malloc(sizeof(double) * n);
  lint N = n, lwork = 2*n + n*n, lrwork = 1 + 5*n + 2*n*n, liwork = 3 + 5*n, info = 0, itype = 1;
  zc *work = malloc(sizeof(zc) * lwork);
  double *rwork = malloc(sizeof(double) * lrwork);
  lint *iwork = malloc(sizeof(lint) * liwork);
  char jobz = 'V';
  const char *name[] = {"zheev ", "zhegv ", "zheevd"};
  build(A, n, band, nzero);
  memcpy(V, A, sizeof(zc) * n * n);
  if (routine == 0)
  {
    zheev_(&jobz, &uplo, &N, V, &N, w, work, &lwork, rwork, &info);
  }
  else if (routine == 1)
  {
    memset(B, 0, sizeof(zc) * n * n);
    for (int i = 0; i < n; i++) RE(B[i + i*n]) = 1;
    zhegv_(&itype, &jobz, &uplo, &N, V, &N, B, &N, w, work, &lwork, rwork, &info);
  }
  else
  {
    zheevd_(&jobz, &uplo, &N, V, &N, w, work, &lwork, rwork, &lrwork, iwork, &liwork, &info);
  }
  double orth, res;
  check(A, V, w, n, &orth, &res);
  printf("%s UPLO=%c  N=%3d  zero trailing rows=%d  INFO=%d   max|V^H V - 1| = %8.1e   max|A V - V W| = %8.1e   %s\n",
         name[routine], uplo, n, nzero, (int)info, orth, res, (orth > 1e-10 || res > 1e-10) ? "<-- WRONG" : "ok");
  free(A); free(V); free(B); free(w); free(work); free(rwork); free(iwork);
}

int main(void)
{
  printf("%s\n\n", interface_name);
  run(0, 'U', 129, 11, 2);   /* N = 129: correct */
  run(0, 'U', 130, 11, 2);   /* N = 130: wrong   */
  run(0, 'U', 130, 11, 1);   /* one zero row: correct */
  run(0, 'U', 138, 11, 5);
  run(0, 'U', 138, 137, 5);  /* dense: wrong as well */
  run(0, 'U', 300, 11, 2);
  run(1, 'U', 138, 11, 5);   /* zhegv, B = I */
  printf("\nthe same matrices with UPLO = 'L' or zheevd:\n");
  run(0, 'L', 138, 11, 5);
  run(1, 'L', 138, 11, 5);
  run(2, 'U', 138, 11, 5);
  return 0;
}