C Program to Perform LU Decomposition of Any Matrix

This is a C Program to perform LU decomposition of a given matrix. Every square matrix A can be decomposed into a product of a lower triangular matrix L and a upper triangular matrix U, as described in LU decomposition.
A = LU

Here is source code of the C Program to Perform LU Decomposition of any Matrix. The C program is successfully compiled and run on a Linux system. The program output is also shown below.

  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4.  
  5. #define foreach(a, b, c) for (int a = b; a < c; a++)
  6. #define for_i foreach(i, 0, n)
  7. #define for_j foreach(j, 0, n)
  8. #define for_k foreach(k, 0, n)
  9. #define for_ij for_i for_j
  10. #define for_ijk for_ij for_k
  11. #define _dim int n
  12. #define _swap(x, y) { typeof(x) tmp = x; x = y; y = tmp; }
  13. #define _sum_k(a, b, c, s) { s = 0; foreach(k, a, b) s+= c; }
  14.  
  15. typedef double **mat;
  16.  
  17. #define _zero(a) mat_zero(a, n)
  18. void mat_zero(mat x, int n) {
  19.     for_ij
  20.             x[i][j] = 0;
  21. }
  22.  
  23. #define _new(a) a = mat_new(n)
  24. mat mat_new(_dim) {
  25.     mat x = malloc(sizeof(double*) * n);
  26.     x[0] = malloc(sizeof(double) * n * n);
  27.  
  28.     for_i
  29.         x[i] = x[0] + n * i;
  30.     _zero(x);
  31.  
  32.     return x;
  33. }
  34.  
  35. #define _copy(a) mat_copy(a, n)
  36. mat mat_copy(void *s, _dim) {
  37.     mat x = mat_new(n);
  38.     for_ij
  39.             x[i][j] = ((double(*)[n]) s)[i][j];
  40.     return x;
  41. }
  42.  
  43. #define _del(x) mat_del(x)
  44. void mat_del(mat x) {
  45.     free(x[0]);
  46.     free(x);
  47. }
  48.  
  49. #define _QUOT(x) #x
  50. #define QUOTE(x) _QUOT(x)
  51. #define _show(a) printf(QUOTE(a)" =");mat_show(a, 0, n)
  52. void mat_show(mat x, char *fmt, _dim) {
  53.     if (!fmt)
  54.         fmt = "%8.4g";
  55.     for_i {
  56.         printf(i ? "      " : " [ ");
  57.         for_j {
  58.             printf(fmt, x[i][j]);
  59.             printf(j < n - 1 ? "  " : i == n - 1 ? " ]\n" : "\n");
  60.         }
  61.     }
  62. }
  63.  
  64. #define _mul(a, b) mat_mul(a, b, n)
  65. mat mat_mul(mat a, mat b, _dim) {
  66.     mat c = _new(c);
  67.     for_ijk
  68.                 c[i][j] += a[i][k] * b[k][j];
  69.     return c;
  70. }
  71.  
  72. #define _pivot(a, b) mat_pivot(a, b, n)
  73. void mat_pivot(mat a, mat p, _dim) {
  74.     for_ij
  75.         {
  76.             p[i][j] = (i == j);
  77.         }
  78.     for_i {
  79.         int max_j = i;
  80.         foreach(j, i, n)
  81.             if (fabs(a[j][i]) > fabs(a[max_j][i]))
  82.                 max_j = j;
  83.  
  84.         if (max_j != i)
  85.             for_k {
  86.                 _swap(p[i][k], p[max_j][k]);
  87.             }
  88.     }
  89. }
  90.  
  91. #define _LU(a, l, u, p) mat_LU(a, l, u, p, n)
  92. void mat_LU(mat A, mat L, mat U, mat P, _dim) {
  93.     _zero(L);
  94.     _zero(U);
  95.     _pivot(A, P);
  96.  
  97.     mat Aprime = _mul(P, A);
  98.  
  99.     for_i {
  100.         L[i][i] = 1;
  101.     }
  102.     for_ij
  103.         {
  104.             double s;
  105.             if (j <= i) {
  106.                 _sum_k(0, j, L[j][k] * U[k][i], s)
  107.                 U[j][i] = Aprime[j][i] - s;
  108.             }
  109.             if (j >= i) {
  110.                 _sum_k(0, i, L[j][k] * U[k][i], s);
  111.                 L[j][i] = (Aprime[j][i] - s) / U[i][i];
  112.             }
  113.         }
  114.  
  115.     _del(Aprime);
  116. }
  117.  
  118. double A3[][3] = { { 1, 3, 5 }, { 2, 4, 7 }, { 1, 1, 0 } };
  119. double A4[][4] = { { 11, 9, 24, 2 }, { 1, 5, 2, 6 }, { 3, 17, 18, 1 }, { 2, 5,
  120.         7, 1 } };
  121.  
  122. int main() {
  123.     int n = 3;
  124.     mat A, L, P, U;
  125.  
  126.     _new(L);
  127.     _new(P);
  128.     _new(U);
  129.     A = _copy(A3);
  130.     _LU(A, L, U, P);
  131.     _show(A);
  132.     _show(L);
  133.     _show(U);
  134.     _show(P);
  135.     _del(A);
  136.     _del(L);
  137.     _del(U);
  138.     _del(P);
  139.  
  140.     printf("\n");
  141.  
  142.     n = 4;
  143.  
  144.     _new(L);
  145.     _new(P);
  146.     _new(U);
  147.     A = _copy(A4);
  148.     _LU(A, L, U, P);
  149.     _show(A);
  150.     _show(L);
  151.     _show(U);
  152.     _show(P);
  153.     _del(A);
  154.     _del(L);
  155.     _del(U);
  156.     _del(P);
  157.  
  158.     return 0;
  159. }

Output:

$ gcc LUDecomposition.cpp
$ ./a.out
 
A
 
1   3   5
2   4   7
1   1   0
 
L
 
1.00000   0.00000   0.00000
0.50000   1.00000   0.00000
0.50000  -1.00000   1.00000
 
U
 
2.00000   4.00000   7.00000
0.00000   1.00000   1.50000
0.00000   0.00000  -2.00000
 
P
 
0   1   0
1   0   0
0   0   1

Sanfoundry Global Education & Learning Series – 1000 C Programs.

advertisement

Here’s the list of Best Books in C Programming, Data Structures and Algorithms.

advertisement
advertisement
Subscribe to our Newsletters (Subject-wise). Participate in the Sanfoundry Certification contest to get free Certificate of Merit. Join our social networks below and stay updated with latest contests, videos, internships and jobs!

Youtube | Telegram | LinkedIn | Instagram | Facebook | Twitter | Pinterest
Manish Bhojasia - Founder & CTO at Sanfoundry
I’m Manish - Founder and CTO at Sanfoundry. I’ve been working in tech for over 25 years, with deep focus on Linux kernel, SAN technologies, Advanced C, Full Stack and Scalable website designs.

You can connect with me on LinkedIn, watch my Youtube Masterclasses, or join my Telegram tech discussions.

If you’re in your 40s–60s and exploring new directions in your career, I also offer mentoring. Learn more here.