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.

Sanfoundry Certification Contest of the Month is Live. 100+ Subjects. Participate Now!
advertisement
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
Manish Bhojasia, a technology veteran with 20+ years @ Cisco & Wipro, is Founder and CTO at Sanfoundry. He lives in Bangalore, and focuses on development of Linux Kernel, SAN Technologies, Advanced C, Data Structures & Alogrithms. Stay connected with him at LinkedIn.

Subscribe to his free Masterclasses at Youtube & technical discussions at Telegram SanfoundryClasses.