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