-3

This is my code for LU Decomposition Crout Method:

function [L, U] = croutreduction(A)
    [row,column]=size(A);
    L=eye(row,column);

    //A = 3x3
    if row==3 then
        U(1,1)=A(1,1); U(1,2)=A(1,2); U(1,3)=A(1,3);
        L(2,1)=A(2,1)/U(1,1); L(3,1)=A(3,1)/U(1,1);

        U(2,2)=A(2,2)-L(2,1)*U(1,2);
        U(2,3)=A(2,3)-L(2,1)*U(1,3);
        L(3,2)=A(3,2)/U(2,2);

        U(3,3)=A(3,3)-L(3,2)*U(2,3);    
    end

    //A = 4x4
    if column==4 then
        U(1,1)=A(1,1); U(1,2)=A(1,2); U(1,3)=A(1,3); U(1,4)=A(1,4);
        L(2,1)=A(2,1)/U(1,1); L(3,1)=A(3,1)/U(1,1); L(4,1)=A(4,1)/U(1,1);

        U(2,2)=A(2,2)-L(2,1)*U(1,2);
        U(2,3)=A(2,3)-L(2,1)*U(1,3);
        U(2,4)=A(2,4)-L(2,1)*U(1,4);
        L(3,2)=(A(3,2)-L(3,1)*U(1,2))/U(2,2);
        L(4,2)=(A(4,2)-L(4,1)*U(1,2))/U(2,2);

        U(3,3)=A(3,3)-(L(3,1)*U(1,3)+L(3,2)*U(2,3));  
        U(3,4)=A(3,4)-(L(3,1)*U(1,4)+L(3,2)*U(2,4));  
        L(4,3)=(A(4,3)-(L(4,1)*U(1,3)+L(4,2)*U(2,3)))/U(3,3);

        U(4,4)=A(4,4)-(L(4,1)*U(1,4)+L(4,2)*U(2,4)+L(4,3)*U(3,4));
    end   
endfunction

How can I modify my code to work with matrices with different dimensions? As you see the code above is only for 3x3 and 4x4 matrix.

Attila
  • 615
  • 4
  • 10
Elra Ghifary
  • 113
  • 1
  • 14
  • Crout reduction code is avaialbe via Google and regarding first question please do attempt to do it first and after you will have any issues and questions we will be more then happy to help you – madbitloman Apr 04 '15 at 21:47
  • Can you check my new code? – Elra Ghifary Apr 06 '15 at 16:33
  • If you have N-dimensional matrix I would suggest you to reshape it into square 2D matrix and perform inversion and then reshape back. – madbitloman Apr 06 '15 at 17:49

1 Answers1

1

You should use for loops instead of hardcoded indices. Based on this example: http://en.wikipedia.org/wiki/Crout_matrix_decomposition I modified the code for Scilab (the original code is for C and Matlab/Octave):

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
endfunction

Hovewer this gives different result than your code, and I didn't check which one is wrong and where...

Attila
  • 615
  • 4
  • 10