Solving for Ax=B where x>=0

[Farfalo]

Hi guys,

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?

Thank you

[fletchjp]

I think the problem is most likely somewhere in your data.

Try checking a hand solution of the problem.

John

[CyLith]

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.

[Julien Langou]

Hello,

Yep, there is no such routine in LAPACK. You want to do something called "nonnegative least squares".

The reference algorithm is given in Lawson, C.L. and R.J. Hanson, Solving Least Squares Problems, Prentice-Hall, 1974, Chapter 23, p. 161.

Here is a more recent code for it: http://hesperia.gsfc.nasa.gov/~schmahl/nnls/nnls.for

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.html
See: Yuancheng Luo and Ramani Duraiswami, "Efficient Parallel Non-Negative Least Squares on Multicore Architectures", SIAM Journal on Scientific Computing 2011.

Cheers,
Julien.

[serkanulu]

This is useful, "nonnegative least squares" (http://hesperia.gsfc.nasa.gov/~schmahl/nnls/nnls.for). i firstly understood, "least squares" and i thinking, imaginary values for example (ax + bi). i^2 (i square) is equal to -1 (minus one). But it's my fault, you defined (x >= 0). i think sleepy.