I have been looking at LAPACK for a few days now, and I have it working to an extent, however I can't seem to find a function that does exactly what I need. I need to solve Ax=B, where I have a 6x6 matrix A with one right-hand side vector and keep x>=0. I was able to use DGESV in order to solve for x, but the simplest solution in most cases involves some negative numbers, which in my case causes some issues as it involves chemical additions to water (you can't add negative amounts! ).

I was wondering if there is a function in LAPACK or a way to set DGESV up differently so that x is greater than or equal to 0?

If your A matrix is square and full rank, then there is only a single unique solution, and if the solution you get from Lapack has negative entries, then as John said, the problem is with your data being inconsistent with expectations. If your A matrix is rank deficient (has a null space), then there are an infinite number of solutions, and you will have to select one by some additional criteria (like positivity). You can find the least norm solution using the built in Lapack routines, and also find the nullspace of A with a QR decomposition; together these two things let you parameterize the entire space of solutions.

And I did not know about this but just found that there is a CUDA / Open MP implementation for the code:http://www.cs.umd.edu/~yluo1/Projects/NNLS.htmlSee: Yuancheng Luo and Ramani Duraiswami, "Efficient Parallel Non-Negative Least Squares on Multicore Architectures", SIAM Journal on Scientific Computing 2011.