alvaroabascar / mad Goto Github PK

B can be just a vector. In the case of gaussj or gaussjbs, adding the capability to solve at the same time for several right hand sides makes sense. In the case of ludcmp we can do the LU decomposition once, and then solve for every right hand side, one at a time. Thus, B should be just a vector.

Current version uses full pivoting but it's not correctly implemented, as it doesn't invert the matrix. Partial pivoting is easier and should work fine, and full pivoting can be fully implemented later, if needed.

It just makes no sense to do it by reference...

In the array "scaling" found in gaussjbs and ludcmp, once pivoting is done in the matrix, the corresponding elements should be interchanged. Otherwise we end up with an array containing false scales.

An alternative strategy to interchanging all the elements in two rows, such as using an indexing array (Numerical Mathematics and Computing By E. Cheney, David Kincaid, section 2.2), or using pointers could improve a lot the efficiency of the algorithms.

Returns wrong answer.

In all the algorithms: gaussj, gaussjbs and ludcmp.