Assert doesn't work with std::numeric limits epsilon, so new eps manu…

[sdharmin]

…ally set

working gauss elimination with LU decomposition

Signed-off-by: Dharmin Shah [email protected]

[JusticeBoi]

Thanks for the epsilon correction, i guess it was ambitious to hope for that precise solution.

For your solution, you can think about a way to solve it without explicitly creating L and U matrices. You can iterate the matrix in a way to only get the upper or lower part. We are looking for a solution with smallest memory footprint.

[sdharmin]

@sariug Yes, I thought it had to be like LU decomposition, but of course I can send the regular Gauss elimation, brute approach

[sariug]

True that I encouraged for LU in the beginning, due to the pseudo code the class was checking as well as an old code I wrote in 2014:

if method  == 1                 % Backslash
    P(:)  = A\s(:);
elseif method == 2
    P(pp) = (Up\(Lp\s(pp)));      % LU
end

The backslash is the solver of MATLAB here. As you can see it first solves for low then up therefore the solver is the same but matrices are changing.

The point I was trying to make there is that the Gauss should solve it no matter what, but the LU decomposition may increase efficiency. Not to write a code depending on LU.

[JusticeBoi]

If the point is to solve the problem, I don't think it is needed to allocate memory for two matrices, it would scale badly. I agree that LU decomposition is also a nice algorithm to be explored. I am just not sure if it is expected behavior of a gauss-elimination based solver to explicitly create L and U matrices. You can gladly prove me wrong by benchmarking it though! Create a new function and test it, might be fun.

[sdharmin]

@sariug added normal Gauss as suggested, hopefully this fulfills the purppose. @JusticeBoi True if the point is solving then no need to create 2 matrices, i think where LU decomposition becomes beneficial is when you have to use the same matrix multiple times with different rhs, imagine a fluid flow with changing dofs and you have to do Ax=b for different b and find the x, at that time the complexity of LU is just n^2 while Gauss remains n^3. Of course for a trivial case like 3*3 it doesn't make sense! Nonetheless, I wanted to give it a shot. Also, I am not sure how to test it...as in the speed, will using multiple rhs give result always or it might fail, I don't know. Let me know about what you think.