/******************************************************************************
gsl_linalg_complex_cholesky_invert()
  Compute the inverse of an Hermitian positive definite matrix in
  Cholesky form.

Inputs: LLT - matrix in cholesky form on input
              A^{-1} = L^{-H} L^{-1} on output

Return: success or error
******************************************************************************/

int
gsl_linalg_complex_cholesky_invert(gsl_matrix_complex * LLT)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else
    {
      size_t N = LLT->size1;
      size_t i, j;
      gsl_vector_complex_view v1;

      /* invert the lower triangle of LLT */
      for (i = 0; i < N; ++i)
        {
          double ajj;
          gsl_complex z;

          j = N - i - 1;

          ajj = 1.0 / GSL_REAL(gsl_matrix_complex_get(LLT, j, j));
          GSL_SET_COMPLEX(&z, ajj, 0.0);
          gsl_matrix_complex_set(LLT, j, j, z);
          ajj = -GSL_REAL(gsl_matrix_complex_get(LLT, j, j));

          if (j < N - 1)
            {
              gsl_matrix_complex_view m;
              
              m = gsl_matrix_complex_submatrix(LLT, j + 1, j + 1,
                                       N - j - 1, N - j - 1);
              v1 = gsl_matrix_complex_subcolumn(LLT, j, j + 1, N - j - 1);

              gsl_blas_ztrmv(CblasLower, CblasNoTrans, CblasNonUnit,
                             &m.matrix, &v1.vector);

              gsl_blas_zdscal(ajj, &v1.vector);
            }
        } /* for (i = 0; i < N; ++i) */

      /*
       * The lower triangle of LLT now contains L^{-1}. Now compute
       * A^{-1} = L^{-H} L^{-1}
       *
       * The (ij) element of A^{-1} is column i of conj(L^{-1}) dotted into
       * column j of L^{-1}
       */

      for (i = 0; i < N; ++i)
        {
          gsl_complex sum;
          for (j = i + 1; j < N; ++j)
            {
              gsl_vector_complex_view v2;
              v1 = gsl_matrix_complex_subcolumn(LLT, i, j, N - j);
              v2 = gsl_matrix_complex_subcolumn(LLT, j, j, N - j);

              /* compute Ainv[i,j] = sum_k{conj(Linv[k,i]) * Linv[k,j]} */
              gsl_blas_zdotc(&v1.vector, &v2.vector, &sum);

              /* store in upper triangle */
              gsl_matrix_complex_set(LLT, i, j, sum);
            }

          /* now compute the diagonal element */
          v1 = gsl_matrix_complex_subcolumn(LLT, i, i, N - i);
          gsl_blas_zdotc(&v1.vector, &v1.vector, &sum);
          gsl_matrix_complex_set(LLT, i, i, sum);
        }

      /* copy the Hermitian upper triangle to the lower triangle */

      for (j = 1; j < N; j++)
        {
          for (i = 0; i < j; i++)
            {
              gsl_complex z = gsl_matrix_complex_get(LLT, i, j);
              gsl_matrix_complex_set(LLT, j, i, gsl_complex_conjugate(z));
            }
        } 

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_invert() */


/*----------------
Testing Functions
----------------*/

int
test_choleskyc_invert_dim(const gsl_matrix_complex * m, double eps)
{
  int s = 0;
  unsigned long i, j, N = m->size1;
  gsl_complex af, bt;
  gsl_matrix_complex * v  = gsl_matrix_complex_alloc(N, N);
  gsl_matrix_complex * c  = gsl_matrix_complex_alloc(N, N);

  gsl_matrix_complex_memcpy(v, m);

  s += gsl_linalg_complex_cholesky_decomp(v);
  s += gsl_linalg_complex_cholesky_invert(v);

  GSL_SET_COMPLEX(&af, 1.0, 0.0);
  GSL_SET_COMPLEX(&bt, 0.0, 0.0);
  gsl_blas_zhemm(CblasLeft, CblasUpper, af, m, v, bt, c);

  /* c should be the identity matrix */

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          int foo;
          gsl_complex cij = gsl_matrix_complex_get(c, i, j);
          double expected, actual;

          /* check real part */
          if (i == j)
            expected = 1.0;
          else
            expected = 0.0;
          actual = GSL_REAL(cij);
          foo = check(actual, expected, eps);
          if (foo)
            printf("REAL(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", 
                   N, N, i,j, actual, expected);
          s += foo;

          /* check imaginary part */
          expected = 0.0;
          actual = GSL_IMAG(cij);
          foo = check(actual, expected, eps);
          if (foo)
            printf("IMAG(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", 
                   N, N, i,j, actual, expected);
          s += foo;

        }
    }

  gsl_matrix_complex_free(v);
  gsl_matrix_complex_free(c);

  return s;
}

int
test_choleskyc_invert(void)
{
  int f;
  int s = 0;

  double dat2[] = { 
      92.303, 0.000,    10.858, 1.798,    
      10.858, -1.798,    89.027, 0.000 
  }; 

  double dat3[] = { 
      59.75,0,       49.25,172.25, 66.75,-162.75,
      49.25,-172.25, 555.5,0,      -429,-333.5,
      66.75,162.75,  -429,333.5,   536.5,0 
  };

  double dat4[] = { 
      102.108, 0.000,    14.721, 1.343,    -17.480, 15.591,    3.308, -2.936,    
      14.721, -1.343,    101.970, 0.000,    11.671, -6.776,    -5.009, -2.665,    
      -17.480, -15.591,    11.671, 6.776,    105.071, 0.000,    3.396, 6.276,    
      3.308, 2.936,    -5.009, 2.665,    3.396, -6.276,    107.128, 0.000 
  }; 

  gsl_matrix_complex_view rv2 = gsl_matrix_complex_view_array(dat2, 2, 2);
  f = test_choleskyc_invert_dim(&rv2.matrix, 2 * 8.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  choleskyc_invert 2x2 Hermitian");
  s += f;

  gsl_matrix_complex_view rv3 = gsl_matrix_complex_view_array(dat3, 3, 3);
  f = test_choleskyc_invert_dim(&rv3.matrix, 2 * 1024.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  choleskyc_invert 3x3 Hermitian");
  s += f;

  gsl_matrix_complex_view rv4 = gsl_matrix_complex_view_array(dat4, 4, 4);
  f = test_choleskyc_invert_dim(&rv4.matrix, 2 * 64.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  choleskyc_invert 4x4 Hermitian");
  s += f;

  return s;
}