Doolittle Algorithm : LU Decomposition

In numerical analysis and linear algebra, LU decomposition (where ‘LU’ stands for ‘lower upper’, and also called LU factorization) factors a matrix as the product of a lower triangular matrix and an upper triangular matrix. Computers usually solve square systems of linear equations using the LU decomposition, and it is also a key step when inverting a matrix, or computing the determinant of a matrix. The LU decomposition was introduced by mathematician Tadeusz Banachiewicz in 1938.

Let A be a square matrix. An LU factorization refers to the factorization of A, with proper row and/or column orderings or permutations, into two factors, a lower triangular matrix L and an upper triangular matrix U, A=LU.

Doolittle Algorithm:

It is always possible to factor a square matrix into a lower triangular matrix and an upper triangular matrix. That is, [A] = [L][U]
Doolittle’s method provides an alternative way to factor A into an LU decomposition without going through the hassle of Gaussian Elimination.
For a general nÃ—n matrix A, we assume that an LU decomposition exists, and write the form of L and U explicitly. We then systematically solve for the entries in L and U from the equations that result from the multiplications necessary for A=LU.

Terms of U matrix are given by:

And the terms for L matrix:

Example

Input :

Output :

Implementation:

C++

 // C++ Program to decompose a matrix into // lower and upper triangular matrix #include  using namespace std;   const int MAX = 100;   void luDecomposition(int mat[][MAX], int n) {     int lower[n][n], upper[n][n];     memset(lower, 0, sizeof(lower));     memset(upper, 0, sizeof(upper));       // Decomposing matrix into Upper and Lower     // triangular matrix     for (int i = 0; i < n; i++)      {         // Upper Triangular         for (int k = i; k < n; k++)         {             // Summation of L(i, j) * U(j, k)             int sum = 0;             for (int j = 0; j < i; j++)                 sum += (lower[i][j] * upper[j][k]);               // Evaluating U(i, k)             upper[i][k] = mat[i][k] - sum;         }           // Lower Triangular         for (int k = i; k < n; k++)          {             if (i == k)                 lower[i][i] = 1; // Diagonal as 1             else             {                 // Summation of L(k, j) * U(j, i)                 int sum = 0;                 for (int j = 0; j < i; j++)                     sum += (lower[k][j] * upper[j][i]);                   // Evaluating L(k, i)                 lower[k][i]                     = (mat[k][i] - sum) / upper[i][i];             }         }     }       // setw is for displaying nicely     cout << setw(6)           << "      Lower Triangular"          << setw(32)          << "Upper Triangular" << endl;       // Displaying the result :     for (int i = 0; i < n; i++)      {         // Lower         for (int j = 0; j < n; j++)             cout << setw(6) << lower[i][j] << "\t";         cout << "\t";           // Upper         for (int j = 0; j < n; j++)             cout << setw(6) << upper[i][j] << "\t";         cout << endl;     } }   // Driver code int main() {     int mat[][MAX]         = { { 2, -1, -2 }, { -4, 6, 3 }, { -4, -2, 8 } };       luDecomposition(mat, 3);     return 0; }

Java

 // Java Program to decompose a matrix into // lower and upper triangular matrix class GFG {       static int MAX = 100;     static String s = "";     static void luDecomposition(int[][] mat, int n)     {         int[][] lower = new int[n][n];          int[][] upper = new int[n][n];           // Decomposing matrix into Upper and Lower         // triangular matrix         for (int i = 0; i < n; i++)         {             // Upper Triangular             for (int k = i; k < n; k++)              {                 // Summation of L(i, j) * U(j, k)                 int sum = 0;                 for (int j = 0; j < i; j++)                     sum += (lower[i][j] * upper[j][k]);                   // Evaluating U(i, k)                 upper[i][k] = mat[i][k] - sum;             }               // Lower Triangular             for (int k = i; k < n; k++)              {                 if (i == k)                     lower[i][i] = 1; // Diagonal as 1                 else                 {                     // Summation of L(k, j) * U(j, i)                     int sum = 0;                     for (int j = 0; j < i; j++)                         sum += (lower[k][j] * upper[j][i]);                       // Evaluating L(k, i)                     lower[k][i]                         = (mat[k][i] - sum) / upper[i][i];                 }             }         }           // setw is for displaying nicely         System.out.println(setw(2) + "     Lower Triangular"                            + setw(10) + "Upper Triangular");           // Displaying the result :         for (int i = 0; i < n; i++)          {             // Lower             for (int j = 0; j < n; j++)                 System.out.print(setw(4) + lower[i][j]                                  + "\t");             System.out.print("\t");               // Upper             for (int j = 0; j < n; j++)                 System.out.print(setw(4) + upper[i][j]                                  + "\t");             System.out.print("\n");         }     }     static String setw(int noOfSpace)     {         s = "";         for (int i = 0; i < noOfSpace; i++)             s += " ";         return s;     }       // Driver code     public static void main(String arr[])     {         int mat[][] = { { 2, -1, -2 },                         { -4, 6, 3 },                         { -4, -2, 8 } };           luDecomposition(mat, 3);     } }   /* This code contributed by PrinciRaj1992 */

Python3

 # Python3 Program to decompose # a matrix into lower and # upper triangular matrix MAX = 100     def luDecomposition(mat, n):       lower = [[0 for x in range(n)]              for y in range(n)]     upper = [[0 for x in range(n)]              for y in range(n)]       # Decomposing matrix into Upper     # and Lower triangular matrix     for i in range(n):           # Upper Triangular         for k in range(i, n):               # Summation of L(i, j) * U(j, k)             sum = 0             for j in range(i):                 sum += (lower[i][j] * upper[j][k])               # Evaluating U(i, k)             upper[i][k] = mat[i][k] - sum           # Lower Triangular         for k in range(i, n):             if (i == k):                 lower[i][i] = 1  # Diagonal as 1             else:                   # Summation of L(k, j) * U(j, i)                 sum = 0                 for j in range(i):                     sum += (lower[k][j] * upper[j][i])                   # Evaluating L(k, i)                 lower[k][i] = int((mat[k][i] - sum) /                                   upper[i][i])       # setw is for displaying nicely     print("Lower Triangular\t\tUpper Triangular")       # Displaying the result :     for i in range(n):           # Lower         for j in range(n):             print(lower[i][j], end="\t")         print("", end="\t")           # Upper         for j in range(n):             print(upper[i][j], end="\t")         print("")     # Driver code mat = [[2, -1, -2],        [-4, 6, 3],        [-4, -2, 8]]   luDecomposition(mat, 3)   # This code is contributed by mits

C#

 // C# Program to decompose a matrix into // lower and upper triangular matrix using System;   class GFG {       static int MAX = 100;     static String s = "";     static void luDecomposition(int[, ] mat, int n)     {         int[, ] lower = new int[n, n];         int[, ] upper = new int[n, n];           // Decomposing matrix into Upper and Lower         // triangular matrix         for (int i = 0; i < n; i++)          {             // Upper Triangular             for (int k = i; k < n; k++)              {                 // Summation of L(i, j) * U(j, k)                 int sum = 0;                 for (int j = 0; j < i; j++)                     sum += (lower[i, j] * upper[j, k]);                   // Evaluating U(i, k)                 upper[i, k] = mat[i, k] - sum;             }               // Lower Triangular             for (int k = i; k < n; k++)              {                 if (i == k)                     lower[i, i] = 1; // Diagonal as 1                 else                 {                     // Summation of L(k, j) * U(j, i)                     int sum = 0;                     for (int j = 0; j < i; j++)                         sum += (lower[k, j] * upper[j, i]);                       // Evaluating L(k, i)                     lower[k, i]                         = (mat[k, i] - sum) / upper[i, i];                 }             }         }           // setw is for displaying nicely         Console.WriteLine(setw(2) + "     Lower Triangular"                           + setw(10) + "Upper Triangular");           // Displaying the result :         for (int i = 0; i < n; i++)          {             // Lower             for (int j = 0; j < n; j++)                 Console.Write(setw(4) + lower[i, j] + "\t");             Console.Write("\t");               // Upper             for (int j = 0; j < n; j++)                 Console.Write(setw(4) + upper[i, j] + "\t");             Console.Write("\n");         }     }     static String setw(int noOfSpace)     {         s = "";         for (int i = 0; i < noOfSpace; i++)             s += " ";         return s;     }       // Driver code     public static void Main(String[] arr)     {         int[, ] mat = { { 2, -1, -2 },                         { -4, 6, 3 },                         { -4, -2, 8 } };           luDecomposition(mat, 3);     } }   // This code is contributed by Princi Singh

PHP

 

Javascript

 

Output

      Lower Triangular                Upper Triangular
1         0         0             2        -1        -2
-2         1         0             0         4        -1
-2        -1         1             0         0         3    

If you like GeeksforGeeks and would like to contribute, you can also write an article using write.geeksforgeeks.org or mail your article to review-team@geeksforgeeks.org. See your article appearing on the GeeksforGeeks main page and help other Geeks.

Feeling lost in the world of random DSA topics, wasting time without progress? It's time for a change! Join our DSA course, where we'll guide you on an exciting journey to master DSA efficiently and on schedule.
Ready to dive in? Explore our Free Demo Content and join our DSA course, trusted by over 100,000 geeks!

Previous
Next