Doolittle Algorithm: LU Decomposition

LU decomposition is a matrix factorization technique that decomposes a matrix into the product of a lower triangular matrix (L) and an upper triangular matrix (U). The Doolittle algorithm is a specific type of LU decomposition where the diagonal elements of the lower triangular matrix (L) are all ones. This article explains the Doolittle algorithm and provides a C implementation.

Fundamental Concepts / Prerequisites

To understand the Doolittle algorithm, you should be familiar with the following concepts:

• Matrices: Basic matrix operations (addition, subtraction, multiplication).
• Triangular Matrices: Understanding of lower triangular and upper triangular matrices. A lower triangular matrix has all elements above the main diagonal equal to zero. An upper triangular matrix has all elements below the main diagonal equal to zero.
• Linear Algebra: Basic principles of linear algebra.

Implementation in C

#include <stdio.h>
#include <stdlib.h>

// Function to perform Doolittle LU decomposition
void luDecomposition(double** A, int n, double** L, double** U) {
    for (int i = 0; i < n; i++) {
        // Upper Triangular Matrix (U)
        for (int k = i; k < n; k++) {
            double sum = 0;
            for (int j = 0; j < i; j++) {
                sum += (L[i][j] * U[j][k]);
            }
            U[i][k] = A[i][k] - sum;
        }

        // Lower Triangular Matrix (L)
        for (int k = i; k < n; k++) {
            if (i == k) {
                L[i][i] = 1; // Diagonal elements of L are 1
            } else {
                double sum = 0;
                for (int j = 0; j < i; j++) {
                    sum += (L[k][j] * U[j][i]);
                }
                L[k][i] = (A[k][i] - sum) / U[i][i];
            }
        }
    }
}

// Function to print a matrix
void printMatrix(double** matrix, int n) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            printf("%lf\t", matrix[i][j]);
        }
        printf("\n");
    }
}

int main() {
    int n = 3;  // Size of the matrix

    // Example matrix A
    double** A = (double**)malloc(n * sizeof(double*));
    for (int i = 0; i < n; i++) {
        A[i] = (double*)malloc(n * sizeof(double));
    }
    A[0][0] = 2; A[0][1] = 1; A[0][2] = 1;
    A[1][0] = 4; A[1][1] = 1; A[1][2] = 0;
    A[2][0] = -2; A[2][1] = 2; A[2][2] = 1;

    // Allocate memory for L and U
    double** L = (double**)malloc(n * sizeof(double*));
    for (int i = 0; i < n; i++) {
        L[i] = (double*)malloc(n * sizeof(double));
    }

    double** U = (double**)malloc(n * sizeof(double*));
    for (int i = 0; i < n; i++) {
        U[i] = (double*)malloc(n * sizeof(double));
    }

    // Perform LU decomposition
    luDecomposition(A, n, L, U);

    // Print L and U matrices
    printf("Lower Triangular Matrix (L):\n");
    printMatrix(L, n);

    printf("\nUpper Triangular Matrix (U):\n");
    printMatrix(U, n);

    // Free allocated memory
    for (int i = 0; i < n; i++) {
        free(A[i]);
        free(L[i]);
        free(U[i]);
    }
    free(A);
    free(L);
    free(U);

    return 0;
}

Code Explanation

The C code implements the Doolittle LU decomposition algorithm. Here's a step-by-step breakdown:

1. Header Files: stdio.h is included for standard input/output operations, and stdlib.h is included for memory allocation.

2. luDecomposition Function:

- Takes the input matrix A, its size n, and pointers to the L and U matrices as arguments.

- The outer loop (i) iterates through rows/columns.

- The first inner loop (k) calculates the elements of the upper triangular matrix U.

- It calculates U[i][k] by subtracting the sum of products of elements from the L and U matrices from A[i][k].

- The second inner loop (k) calculates the elements of the lower triangular matrix L.

- If i == k (diagonal element), L[i][i] is set to 1 (Doolittle's condition).

- Otherwise, it calculates L[k][i] using a similar formula, involving the sum of products from L and U, and dividing by the corresponding element of U (U[i][i]).

3. printMatrix Function: A utility function to print a matrix to the console for verification.

4. main Function:

- Defines the size of the matrix (n = 3).

- Creates a sample matrix A and allocates memory for it and for the L and U matrices.

- Calls the luDecomposition function to perform the decomposition.

- Calls the printMatrix function to display the resulting L and U matrices.

- Frees the dynamically allocated memory to prevent memory leaks.

Complexity Analysis

Time Complexity: The Doolittle algorithm has a time complexity of O(n³), where n is the size of the matrix. This is because the algorithm involves three nested loops. The outer loop iterates `n` times, and the two inner loops also iterate up to `n` times. The dominant operation is the calculation of the sum of products within the inner loops, which takes O(n) time. Therefore, the total time complexity is approximately n * n * n = O(n³).

Space Complexity: The algorithm requires O(n²) space to store the original matrix A, the lower triangular matrix L, and the upper triangular matrix U. Each of these matrices has n² elements. Therefore, the space complexity is O(n²).

Alternative Approaches

Crout Algorithm: The Crout algorithm is another variant of LU decomposition where the diagonal elements of the upper triangular matrix (U) are all ones. This is an alternative to Doolittle's algorithm, which sets the diagonal of the lower triangular matrix (L) to ones. The choice between Crout and Doolittle often depends on specific application requirements. Crout's algorithm also has a time complexity of O(n³) and a space complexity of O(n²), similar to Doolittle's algorithm.

Conclusion

The Doolittle algorithm provides a method to decompose a matrix into lower and upper triangular matrices, which can be used to solve linear systems of equations efficiently. This article provided a C implementation of the Doolittle algorithm, along with explanations of the code and complexity analysis. Understanding LU decomposition and its variants, such as the Doolittle algorithm, is crucial in various fields, including numerical analysis and scientific computing.