/**
 * Solving Ax = b using Cholesky Decomposition in C using LAPACK
 * Compile with: gcc test_chol.c -lblas -llapack -llapacke -o test_chol
 */
#include <stdio.h>
#include <lapacke.h>

int main() {
  // Define matrix A and vector b in A * x = b
  // clang-format off
  double A[9] = {
    2.0, -1.0, 0.0,
    -1.0, 2.0, -1.0,
    0.0, -1.0, 1.0
  };
  double b[3] = {1.0, 0.0, 0.0};
  int m = 3;
  // clang-format on

  // Factor matrix A using Cholesky decomposition
  int info = 0;
  int lda = m;
  int n = m;
  char uplo = 'U';
  // dpotrf_(&uplo, &n, a, &lda, &info);
  info = LAPACKE_dpotrf(LAPACK_ROW_MAJOR, uplo, n, A, lda);
  if (info != 0) {
    printf("Failed to decompose A using Cholesky Decomposition!\n");
  }

  // Solve Ax = b using Cholesky decomposed A from above
  // Note: As per documentation the solution is written to vector b
  int nhrs = 1;
  int ldb = m;
  // dpotrs_(&uplo, &n, &nhrs, a, &lda, x, &ldb, &info);
  info = LAPACKE_dpotrs(LAPACK_ROW_MAJOR, uplo, n, nhrs, A, lda, b, ldb);
  if (info != 0) {
    printf("Failed to solve Ax = b!\n");
  }

  // Print solution: x should be [1; 1; 1]
  printf("x: [%f, %f, %f]\n", b[0], b[1], b[2]);

  return 0;
}