Skip to content

Instantly share code, notes, and snippets.

@chutsu
Last active March 16, 2023 12:38
Show Gist options
  • Save chutsu/7ba29d2e9b7150ccbf1e241e33a2e310 to your computer and use it in GitHub Desktop.
Save chutsu/7ba29d2e9b7150ccbf1e241e33a2e310 to your computer and use it in GitHub Desktop.
Solving Ax = b using Cholesky Decomposition in C using LAPACK
/**
* 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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment