Last active
March 16, 2023 12:38
-
-
Save chutsu/7ba29d2e9b7150ccbf1e241e33a2e310 to your computer and use it in GitHub Desktop.
Solving Ax = b using Cholesky Decomposition in C using LAPACK
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/** | |
* 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