Crout matrix decomposition

In linear algebra, the Crout matrix decomposition is an LU decomposition which decomposes a matrix into a lower triangular matrix (L), an upper triangular matrix (U) and, although not always needed, a permutation matrix (P). It was developed by Prescott Durand Crout. [1]

The Crout matrix decomposition algorithm differs slightly from the Doolittle method. Doolittle's method returns a unit lower triangular matrix and an upper triangular matrix, while the Crout method returns a lower triangular matrix and a unit upper triangular matrix.

So, if a matrix decomposition of a matrix A is such that:

A = LDU

being L a unit lower triangular matrix, D a diagonal matrix and U a unit upper triangular matrix, then Doolittle's method produces

A = L(DU)

and Crout's method produces

A = (LD)U.

Implementations

C implementation:

void crout(double const **A, double **L, double **U, int n) {
	int i, j, k;
	double sum = 0;

	for (i = 0; i < n; i++) {
		U[i][i] = 1;
	}

	for (j = 0; j < n; j++) {
		for (i = j; i < n; i++) {
			sum = 0;
			for (k = 0; k < j; k++) {
				sum = sum + L[i][k] * U[k][j];	
			}
			L[i][j] = A[i][j] - sum;
		}

		for (i = j; i < n; i++) {
			sum = 0;
			for(k = 0; k < j; k++) {
				sum = sum + L[j][k] * U[k][i];
			}
			if (L[j][j] == 0) {
				printf("det(L) close to 0!\n Can't divide by 0...\n");
				exit(EXIT_FAILURE);
			}
			U[j][i] = (A[j][i] - sum) / L[j][j];
		}
	}
}

Octave/Matlab implementation:

   function [L, U] = LUdecompCrout(A)
        
        [R, C] = size(A);
        for i = 1:R
            L(i, 1) = A(i, 1);
            U(i, i) = 1;
        end
        for j = 2:R
            U(1, j) = A(1, j) / L(1, 1);
        end
        for i = 2:R
            for j = 2:i
                L(i, j) = A(i, j) - L(i, 1:j - 1) * U(1:j - 1, j);
            end
            
            for j = i + 1:R
                U(i, j) = (A(i, j) - L(i, 1:i - 1) * U(1:i - 1, j)) / L(i, i);
            end
        end
   end
gollark: It's because your programs end with `error`, presumably.
gollark: Hmm, so what you're saying is that tail calls are incomprehensible dark magic.
gollark: Tail calls are where you call a function at the end of a function or something, and this is magically optimized better because something something stack.
gollark: Couldn't you just PR it to not not do that?
gollark: You can, as far as I know, emulate pcall-type stuff with temporary coroutines (which is very hacky but oh well), and those would probably not be subject to stack stuff.

References
  1. Press, William H. (2007). Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge University Press. pp. 50–52. ISBN 9780521880688.
This article is issued from Wikipedia. The text is licensed under Creative Commons - Attribution - Sharealike. Additional terms may apply for the media files.