Simplex routine in NR in C, pp.440


aliboray
03-11-2002, 10:06 PM
Sorry in case this is reposted as I am new to this forum.

I have the following question: why is it that the first for loop below on pp.440 starts from m1+1 instead of from 1? After all, starting from m1+1 would only compute the auxiliary objective function for the last two rows when m1=2, m=4 as per the example in the NR in C example book on the Simplex method???

if (m2+m3) {
for (i=m1+1;i<=m2;i++) l3[i]=1;
for (k=1;k<=(n+1);k++)
{
q1=0.0;
for (i=1;i<=m;i++)
{
q1 += a[i+1][k];
}
a[m+2][k] = -q1;
}
:confused:

Bill Press
03-11-2002, 10:42 PM
This is a bug that was already corrected in a previous reprinting. You are correct that the lower limit should be 1.

I think there were also other changes in this routine. You can view updated code (version 2.08) in the Books-On-Line resource:

http://lib-www.lanl.gov/numerical/bookcpdf/c10-8.pdf

Thanks for your post!

aliboray
03-11-2002, 11:22 PM
Thanks for your post Prof. Press,

I notice in the post that the revised code runs as follows:

if (m2+m3) {
for (i=1;i<=m2;i++) l3[i]=1;
for (k=1;k<=(n+1);k++)
{
q1=0.0;
for (i=m1+1;i<=m;i++)
{
q1 += a[i+1][k];
}
a[m+2][k] = -q1;
}


For the part below, I am still under some confusion. In the text to the section, it is mentioned that each entry in the terminal row,
a[m+2][j=2...n+1] is the SUM of all column entries above it, not just column entries above it starting from row m1+1 to m.


...
for (i=m1+1;i<=m;i++)
{
q1 += a[i+1][k];
}
....


I checked your reference on Kunzi adn his Algol code confirms this. Is it that the auxilliary row is NOT as given by equation (10.8.17)???

Any clarification would be very welcome.
Thanks.

Saul Teukolsky
03-12-2002, 09:26 AM
Yes, the text does not correspond exactly to the code that is given. (Not just at the point you mention, but in several places.) The text is meant to be basically pedagogical, and is fairly close to the code.

The basic point here is that "less than" constraints are easily handled to provide an initial feasible basic vector. You don't really need the machinery of introducing artificial variables. You only need the slack variables. You get a feasible basic vector by setting all the *original* variables in those constraints to zero. You only need the "phase one" step to handle the "greater than" and the "equals" constraints. Accordingly, you can construct the objective function for phase one by omitting the first m1 constraints.

There are many other variations on how to do phase one - the procedure in our code is by no means unique.

Let me add that our simplx code is rather dated. More modern approaches formulate the problem in terms of matrices. You can read about this in some of the textbooks we reference.

aliboray
03-14-2002, 01:07 AM
Hi,

Thank you very much for your earlier responses.
This is regarding the simplex code on pp.441 of NR in C.

Am I correct in assuming that for the example given (system 10.8.15), the initial basis is (y1, y2, -y3, z4), since you say on pp.440 as a comment to the code:

"...m2 type constraints have their slack variable initially left hand, with a minus sign, and their artificial variable handled implicitly during the first exhange. m3 type constraints have their artificial variable initially left hand".

Then I can undersand that you multiply the pivot column (if the pivot element corresponds to an m2 type contraint) by -1 in the following section of code on pp441:


else
{
kh=iposv[ip]-m1-n;
if (kh >= 1 && l3[kh])
{
l3[kh]=0;
++a[m+2][kp+1];
for (i=1;i<=m+2;i++)
a[i][kp+1] = -a[i][kp+1];

}
}

This is due to the fact that -y3, in the first exhange, exits the basis, and the tableau is now written in terms of the non-basis variables (x1,x2,x3, -y3) on the right hand side, instead of (x1, x2, x3, +y3), am I correct??

What I don't understand is why you increment the objective function entry corresponding to the pivot element, i.e. the code fragment above:
++a[m+2][kp+1];

I admire the virtuosity of your implementation, but am running into conceptual difficulties which I can't find answers to in the references you cite, hence an answer would be much appreciated.
Incidentally, might you know if any code is available on the Wolfe's algorithm, relying on Simplex?

Regards,

Ali Bora

Saul Teukolsky
03-14-2002, 08:46 AM
First, you can think of the basis as (y1,y2,y3,z4). I'm not sure what
it means to put a minus sign in front of y3 - you're just making a list of
the variables. And yes, that's why there's the -1 you mentioned. The
value of ther objective function must be updated every time you exchange
variables in the basis. Since there was a minus 1 in the objective
function corresponding to y3, you increment by 1 when y3 exits the basis.

As for Wolfe's algorithm, sorry I don't know of any code.