## digitalmars.D.learn - Linear system solver in D?

- BCS (15/15) Feb 18 2008 I am going to have a system of equations like this
- Christopher Wright (10/32) Feb 18 2008 You definitely want bindings rather than native D. D just hasn't been
- Bill Baxter (8/43) Feb 18 2008 Multiarray has bindings to LAPACK.
- BCS (4/56) Feb 19 2008 thanks, both or you, I'll look at those. Performance isn't that big an
- Bill Baxter (23/81) Feb 19 2008 Are the equations the same each time? I.e. do the a_m1 coeffs stay the
- BCS (9/33) Feb 19 2008 no such luck. This will be part of a newton Newton–Raphson non-linear ...
- Guillaume B. (4/62) Feb 19 2008 lp_solve is nice too: http://lpsolve.sourceforge.net/5.5/ ... and it ca...
- Bill Baxter (14/75) Feb 19 2008 Hmm, I was assuming he was after D code. If not, for instance if you
- Guillaume B. (4/83) Feb 19 2008 Well, lp_solve is in C so the bindings should be fairly easy to write.....
- BCS (4/9) Feb 19 2008 how is it's API? What I'll have is an array of real's in memeory.
- Bill Baxter (12/26) Feb 19 2008 Guiallumes link looks to be for integer programming, not floating point
- Guillaume B. (23/53) Feb 20 2008 It actually supports floating points number... And if your C compiler

I am going to have a system of equations like this a_11*x_1 + a_12*x2 + ... a_1n*x_n = y_1 a_21*x_1 + a_22*x2 + ... a_2n*x_n = y_2 . . . a_m1*x_1 + a_m2*x2 + ... a_mn*x_n = y_m y_* and a_* known, I need to find x_* What is the best available solver for such a system that works under D? C bindings would work, D code would be better and I'd rather stay away from non portable (uses __ asm and has no port to other system). If no one knows of a good lib that is ready to use, what is a good C lib to do bindings for? p.s. I'm going to be putting this in a non-linear root finder, has someone already written on of those for D.

Feb 18 2008

BCS wrote:I am going to have a system of equations like this a_11*x_1 + a_12*x2 + ... a_1n*x_n = y_1 a_21*x_1 + a_22*x2 + ... a_2n*x_n = y_2 .. .. .. a_m1*x_1 + a_m2*x2 + ... a_mn*x_n = y_m y_* and a_* known, I need to find x_* What is the best available solver for such a system that works under D? C bindings would work, D code would be better and I'd rather stay away from non portable (uses __ asm and has no port to other system). If no one knows of a good lib that is ready to use, what is a good C lib to do bindings for? p.s. I'm going to be putting this in a non-linear root finder, has someone already written on of those for D.You definitely want bindings rather than native D. D just hasn't been around long enough for people to make decent math libraries for it; most of the people with the required skills are still transitioning from Fortran. You could use GLPK -- it's a linear solver that accepts a superset of AMPL. If you're doing serious work on large data sets, go with CPLEX. If you manage to write something that does any better than GLPK, start a company. CPLEX is significantly better, but you might be able to make some money if you marketed it toward smaller research projects for $500 or so.

Feb 18 2008

Christopher Wright wrote:BCS wrote:Multiarray has bindings to LAPACK. http://www.dsource.org/projects/multiarray There are bindings for GSL which I think uses LAPACK also somewhere on dsource (either its own project or maybe it was in the 'bindings' project). I'm working on a new d library that will wrap LAPACK and some sparse lib like SuperLU and/or TAUCS. The new lib is based loosely on FLENS. --bbI am going to have a system of equations like this a_11*x_1 + a_12*x2 + ... a_1n*x_n = y_1 a_21*x_1 + a_22*x2 + ... a_2n*x_n = y_2 .. .. .. a_m1*x_1 + a_m2*x2 + ... a_mn*x_n = y_m y_* and a_* known, I need to find x_* What is the best available solver for such a system that works under D? C bindings would work, D code would be better and I'd rather stay away from non portable (uses __ asm and has no port to other system). If no one knows of a good lib that is ready to use, what is a good C lib to do bindings for? p.s. I'm going to be putting this in a non-linear root finder, has someone already written on of those for D.You definitely want bindings rather than native D. D just hasn't been around long enough for people to make decent math libraries for it; most of the people with the required skills are still transitioning from Fortran. You could use GLPK -- it's a linear solver that accepts a superset of AMPL. If you're doing serious work on large data sets, go with CPLEX. If you manage to write something that does any better than GLPK, start a company. CPLEX is significantly better, but you might be able to make some money if you marketed it toward smaller research projects for $500 or so.

Feb 18 2008

Bill Baxter wrote:Christopher Wright wrote:thanks, both or you, I'll look at those. Performance isn't that big an issue as I'm only looking at about 15-30 equations and a few minutes run time would be ok, but I'm going to have to run ~1500 passes through it.BCS wrote:Multiarray has bindings to LAPACK. http://www.dsource.org/projects/multiarray There are bindings for GSL which I think uses LAPACK also somewhere on dsource (either its own project or maybe it was in the 'bindings' project). I'm working on a new d library that will wrap LAPACK and some sparse lib like SuperLU and/or TAUCS. The new lib is based loosely on FLENS. --bbI am going to have a system of equations like this a_11*x_1 + a_12*x2 + ... a_1n*x_n = y_1 a_21*x_1 + a_22*x2 + ... a_2n*x_n = y_2 .. .. .. a_m1*x_1 + a_m2*x2 + ... a_mn*x_n = y_m y_* and a_* known, I need to find x_* What is the best available solver for such a system that works under D? C bindings would work, D code would be better and I'd rather stay away from non portable (uses __ asm and has no port to other system). If no one knows of a good lib that is ready to use, what is a good C lib to do bindings for? p.s. I'm going to be putting this in a non-linear root finder, has someone already written on of those for D.You definitely want bindings rather than native D. D just hasn't been around long enough for people to make decent math libraries for it; most of the people with the required skills are still transitioning from Fortran. You could use GLPK -- it's a linear solver that accepts a superset of AMPL. If you're doing serious work on large data sets, go with CPLEX. If you manage to write something that does any better than GLPK, start a company. CPLEX is significantly better, but you might be able to make some money if you marketed it toward smaller research projects for $500 or so.

Feb 19 2008

BCS wrote:Bill Baxter wrote:Are the equations the same each time? I.e. do the a_m1 coeffs stay the same for all 1500 passes? or do they change each time? If they stay the same then you can factorize A once and do the much faster back-substitution 1500 times. If you don't care about performance just search the web for some C/C++/C#/Java code for Gaussian Elimination with Partial (or better, Full) pivoting. It should be pretty straightforward to port. Then you won't have to worry about dependencies and building BLAS/LAPACK etc. Actually you didn't say, but are m and n not the same? If not then you need to use least-squares. If you have more equations than unknowns then you cannot generally find an exact solution, but you can find the solution that minimizes the 2-norm of the residual ||Ax-b||. If you have too few equations then there are many exact solutions, so generally you try to find one that has the smallest 2-norm. I.e. ||x|| is smallest among all possible solutions. One way to solve a least squares problem is to use the normal equations -- you multiply both sides of the equation by A transpose (A^T): A^T A x = A^T b (A^T A) is square, so you can use a regular linear solver on it (like gaussian elimination). --bbChristopher Wright wrote:thanks, both or you, I'll look at those. Performance isn't that big an issue as I'm only looking at about 15-30 equations and a few minutes run time would be ok, but I'm going to have to run ~1500 passes through it.

Feb 19 2008

Reply to Bill,BCS wrote:no such luck. This will be part of a newton Newton–Raphson non-linear solver for doing a numerical approximation of a system of non-linear differential equations. (a.k.a. bad, worse, nasty)thanks, both or you, I'll look at those. Performance isn't that big an issue as I'm only looking at about 15-30 equations and a few minutes run time would be ok, but I'm going to have to run ~1500 passes through it.Are the equations the same each time? I.e. do the a_m1 coeffs stay the same for all 1500 passes? or do they change each time? If they stay the same then you can factorize A once and do the much faster back-substitution 1500 times.If you don't care about performance just search the web for some C/C++/C#/Java code for Gaussian Elimination with Partial (or better, Full) pivoting. It should be pretty straightforward to port.performance takes a back seat to complexity but is still an issue. OTOH generating bindings for canned C code is less complex than porting it.Then you won't have to worry about dependencies and building BLAS/LAPACK etc.building and what not isn't a big issue.Actually you didn't say, but are m and n not the same?On the way out the door after posting that I relized that I put in NxM rather than NxN. oops.--bb

Feb 19 2008

BCS wrote:Bill Baxter wrote:lp_solve is nice too: http://lpsolve.sourceforge.net/5.5/ ... and it can read many file formats... Guillaume

Feb 19 2008

Guillaume B. wrote:BCS wrote:Hmm, I was assuming he was after D code. If not, for instance if you just have some equations you want to solve, then I'd just write a short python/numpy script: --- from numpy import * A = asarray([[a_11, a_12 ...], [a_21, a_22 ...], ... [a_m1, a_m2 ...]]) y = asarray([y_1, y_2...]) x = lstsq(A,y)[0] --- Voila

Feb 19 2008

Bill Baxter wrote:Guillaume B. wrote:Well, lp_solve is in C so the bindings should be fairly easy to write... I suggested that since others suggested using bindings to existing libraries. GuillaumeBCS wrote:Hmm, I was assuming he was after D code. If not, for instance if you just have some equations you want to solve, then I'd just write a short python/numpy script: --- from numpy import * A = asarray([[a_11, a_12 ...], [a_21, a_22 ...], ... [a_m1, a_m2 ...]]) y = asarray([y_1, y_2...]) x = lstsq(A,y)[0] --- Voila

Feb 19 2008

Reply to Guillaume B.,lp_solve is nice too: http://lpsolve.sourceforge.net/5.5/ ... and it can read many file formats... Guillaumehow is it's API? What I'll have is an array of real's in memeory. Hmm that's a though, does anyone know if any of the listed libs work with 80bit FP?

Feb 19 2008

BCS wrote:Reply to Guillaume B.,Guiallumes link looks to be for integer programming, not floating point calcs. LAPACK or something that wraps it is what you want. I highly recommend you look at the BLASLAPACK dir of dsource/projects/multiarray. The LAPACK wrappers have been tested on both Windows and Ubuntu. I have instructions in both dirs for how to install in either case, and there's a simple test program there that does a linear solve on a small matrix. http://www.dsource.org/projects/multiarray/browser/trunk/multiarray/BLASLAPACK/lapacktest.d The blas and lapack are separate from the rest of multiarray, so you can just check out that BLASLAPACK dir and be good to go. --bblp_solve is nice too: http://lpsolve.sourceforge.net/5.5/ ... and it can read many file formats... Guillaumehow is it's API? What I'll have is an array of real's in memeory. Hmm that's a though, does anyone know if any of the listed libs work with 80bit FP?

Feb 19 2008

Bill Baxter wrote:BCS wrote:http://www.dsource.org/projects/multiarray/browser/trunk/multiarray/BLASLAPACK/lapacktest.dReply to Guillaume B.,Guiallumes link looks to be for integer programming, not floating point calcs. LAPACK or something that wraps it is what you want. I highly recommend you look at the BLASLAPACK dir of dsource/projects/multiarray. The LAPACK wrappers have been tested on both Windows and Ubuntu. I have instructions in both dirs for how to install in either case, and there's a simple test program there that does a linear solve on a small matrix.lp_solve is nice too: http://lpsolve.sourceforge.net/5.5/ ... and it can read many file formats... Guillaumehow is it's API? What I'll have is an array of real's in memeory. Hmm that's a though, does anyone know if any of the listed libs work with 80bit FP?The blas and lapack are separate from the rest of multiarray, so you can just check out that BLASLAPACK dir and be good to go. --bbIt actually supports floating points number... And if your C compiler supports 80 bits floating point, you can probably compile it to support it since the functions are like this: unsigned char add_constraint(lprec *lp, REAL *row, int constr_type, REAL rh); Although by default, REAL is probably a typedef for double. And lp_solve is actually a solver for Linear Programs... It's an optimizer, not a simple solver for a system of linear equations... It's to solve things like this: maximize 143x + 60y where 120x + 210y <= 15000 110x + 30y <= 4000 x + y <= 75 x >= 0 y >= 0 So it might not be exactly what you want... You can see an example problem (the explanations are very good) with an example of the C API there: http://lpsolve.sourceforge.net/5.5/formulate.htm Hope this helps! Guillaume

Feb 20 2008