Having problems solving Ax = b using LAPACKE in C

81 views Asked by At

I have an Ax = b problem that I'd like to solve using Cholesky Decomposition, where:

A = [ 2.0, -1.0, 0.0;
     -1.0, 2.0, -1.0;
      0.0, -1.0, 1.0];

b = [1.0; 0.0; 0.0];

Solution x should be:

x = [1; 1; 1];

But I'm having difficulty obtaining the answer using LAPACKE (C interface to LAPACK) and in C. Here is my minimal example code:

// test_chol.c
#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';
  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;
  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;
}

Compile and run with:

gcc test_chol.c -lblas -llapack -llapacke -o test_chol
0

There are 0 answers