I am looking for a FORTRAN routine that can solve a linear system of equations. I don't care if the algorithm is fast or slow, as long as it is "robust", I am happy.

I have tried some of the methods given in Recipes and NetLib including Gaussian ellimination, but I keep getting singularity error, but when I solve the same system of equations in MATLAB using the simple operaor \ (ie, X = A \ B) it solves it and gives the correct answer.

I have spent enough time on this, so I reaaly need to find a nice code so that I can use it and carry on my work.

Thanks very much,

Nora

Check the condition number of the coefficient matrix with MATLAB; also look at the determinant. It could be that the matrix is so poorly conditioned that the Fortran routines classify it as singular, while MATLAB or Mathematica may switch to carrying extra precision or to using a special method to handle ill-conditioned matrices (a "robust" method). Use Fortran routines that use full-pivoting or other reordering technique to minimize the round-off build-up, since you do not mind sacrificing efficiency. Also, be sure that you are supplying the coefficient matrix and the right hand side vector in the correct format, for each routine.

Thanks for the advice. Do you know any links or resources that I can find such a robust routine?

Nora

Nora,

singulary is probably due to having small numbers or zeros on diagonal entries of system, so you have to use pivoting or using methods that handles indefinite systems,

BILUM package of Y. Saad has tihs capability but simpler solution is usnig (i recommend this for u) MA27 routine of HSL (fortran 77, freely availible) for solution of summetric indefinite systems (is sparce direct), or MA47 for non-symmertic indefinite system.

Generally for ill-posed systems sparse direct solvers are recommended, other famous codes are umfpack and PARDISO (for parallel run)