The following program can be used to analyze two player constant sum games. It handles games with multiple states and with some types of recursion. I also have a downloadable version that does not have specials (<, > and &) escaped.
/*
* Game Strategy Solver
*
* Written by Bruno Wolff III
* December 10, 1989
*
* Changed "=-" to "= -" 12/26/91
*
* Made test for suitable pivots slightly larger than 0 rather than against
* exactly 0. 6/4/96
*
* This program will for solve for the value of each state of a game
* and for the optimal strategy for each state. Each state may have
* its own set of strategies and state transitions.
*
* The input to the program is:
* The type of recursion allowed.
* (1 for no recursion, 2 for variable recursion, 3 for total recursion)
* The input format (1 for dense, 2 for sparse).
* The scale value to divide input values by.
* The number of states.
* If recursion allowed:
* The error limit to be used to check for completion.
* The maximum number of iterations.
* For each state:
* The number of row and number of column strategies.
* The immediate payoff table from the row point of view.
* The transition table to allowed states.
*
* The output from the program is:
* If recursion allowed:
* The number of iterations used.
* The maximum change in a state value in the last iteration.
* The value of each state.
* For total recursion:
* The value per turn of the game.
* The probability of being in each state.
* For each state:
* The optimal row strategy mix.
* The optimal column strategy mix.
*
* The problem is solved by treating each state as a seperate game
* using the previous state values when accounting for state transitions.
* Each game is then interpretted as a linear programming problem and
* solved using the simplex method. The new state values are checked against
* the old state values and the number of iterations is checked against the
* maximum number of iterations to see if any more iterations are needed.
* If more iterations are needed then new state values are calculated
* using a multinomial Newton's method (unless I-F' is singular, in which
* case the current new values are used) using the values of the transition
* matrices (evaluted using the latest strategy mixes) as the entries of F'.
* The new state values are then used as the old state values in the next
* iteration. For total recursion an extra condition is added to I-F' so
* that it will not be singular. The condition is that the value of each
* state times the probability of being in that state sum to zero. In this
* mode only one infinite loop of states is allowed.
*/
/*
* Include files.
*/
#include <stdio.h>
/*
* Declare function types.
*/
void main();
void gamsol();
int lineqs();
void badinp();
/*
* Main program.
*/
void main()
{
int rec; /* type of recursion allowed */
int inp; /* format of input */
float scale; /* scale factor for input */
int nstate; /* number of states */
int miter; /* maximum number of iterations */
int iter; /* number of iterations */
float merr; /* maximum error */
float err; /* largest state value change in last iteration */
int *m; /* table of number of row startegies */
int *n; /* table of number of column strategies */
struct entry { /* sparse table entry */
int estate; /* transition to state, 0 for payoff */
float eval; /* value of entry */
struct entry *enext; /* next entry in this state, row, column */
};
float ****dp; /* [state from][state to][row][col] dense payoff table */
struct entry ****sp; /* [state from][row][col] sparse payoff table */
struct entry *ent; /* scratch */
float *old; /* table of old state values */
float *val; /* table of current state values */
float **row; /* [state][strategy] row strategy mixes */
float **col; /* [state][strategy] column strategy mixes */
float **a; /* [row][col] scratch matrix used for gamsol and lineqs calls */
float *b; /* scratch vector used for lineqs call */
float *x; /* scratch vector used for lineqs call */
float v; /* scratch */
int kk; /* used to keep track of allowed state transitions */
int i,j,k,ii,jj; /* scratch */
fprintf(stderr,"\nWhat type of recursion is allowed (1 none, 2 variable, 3 total)?\n");
if (scanf("%d",&rec)<1) badinp();
if (rec<1||rec>3) {
fprintf(stderr,"\nThe type of recursion allowed must be 1, 2 or 3.\n");
exit();
}
fprintf(stderr,"\nWhat input format will be used (1 dense, 2 sparse)?\n");
if (scanf("%d",&inp)<1) badinp();
if (inp<1||inp>2) {
fprintf(stderr,"\nThe input format must be 1 or 2.\n");
exit();
}
fprintf(stderr,"\nWhat is the scale factor for the input?\n");
if (scanf("%f",&scale)<1) badinp();
if (scale==0.) {
fprintf(stderr,"\nThe scale factor must be nonzero.\n");
exit();
}
fprintf(stderr,"\nHow many states does the game have?\n");
if (scanf("%d",&nstate)<1) badinp();
if (nstate < 1) {
fprintf(stderr,"\nThe number of states can not be less than one.\n");
exit();
}
kk=nstate; /* Assume any state transitons are allowed, change later if not */
if (rec==1) {
miter=1; /* Only one iteration needed if there is no recursion */
merr=0.; /* This doesn't matter if there is no recursion */
} else {
fprintf(stderr,"\nWhat is the maximum error?\n");
if (scanf("%f",&merr)<1) badinp();
if (merr<=0.) {
fprintf(stderr,"\nThe maximum error must be greater than zero.\n");
exit();
}
fprintf(stderr,"\nWhat is the maximum number of iterations?\n");
if (scanf("%d",&miter)<1) badinp();
if (miter<1) {
fprintf(stderr,"\nThe maximum number of iterations must be at least one.\n");
exit();
}
}
/*
* Allocate some space for tables.
*/
m=(int *)malloc(sizeof(*m)*nstate); /* number of row strategies table */
n=(int *)malloc(sizeof(*n)*nstate); /* number of column strategies table */
old=(float *)malloc(sizeof(*old)*(nstate+2)); /* old value table */
val=(float *)malloc(sizeof(*val)*(nstate+2)); /* new value table */
row=(float **)malloc(sizeof(*row)*nstate); /* row strategy mixes */
col=(float **)malloc(sizeof(*col)*nstate); /* column strategy mixes */
b=(float *)malloc(sizeof(*b)*(nstate+1));
x=(float *)malloc(sizeof(*x)*(nstate+1));
if (inp==1) { /* Dense input space allocation */
dp=(float ****)malloc(sizeof(*dp)*nstate);
for (i=0; i<nstate; i++) {
if (rec==1) kk=i; /* In case not all transitions are allowed */
dp[i]=(float ***)malloc(sizeof(**dp)*(kk+1));
}
} else { /* Sparse input allocation */
sp=(struct entry ****)malloc(sizeof(*sp)*nstate);
}
/*
* Get number of row and column strategies for each state.
*/
for (i=0; i<nstate; i++) {
fprintf(stderr,"\nWhat is the number of row and the number of column strategies for state %d?\n",i+1);
if (scanf("%d%d",m+i,n+i)<2) badinp();
if (m[i]<1) {
fprintf(stderr,"\nThe number of row strategies must be at least one.\n");
exit();
}
if (n[i]<1) {
fprintf(stderr,"\nThe number of column strategies must be at least one.\n");
exit();
}
row[i]=(float *)malloc(sizeof(**row)*m[i]); /* Space for row strategies */
col[i]=(float *)malloc(sizeof(**col)*n[i]); /* Space for column strategies */
}
/*
* Enter tables using dense input format.
*/
if (inp==1) {
for (i=0; i<nstate; i++) {
if (rec==1) kk=i; /* In case not all transitions are allowed */
for (ii=0; ii<=kk; ii++) { /* Allocate space for each state */
dp[i][ii]=(float **)malloc(sizeof(***dp)*m[i]);
for (k=0; k<m[i]; k++) {
dp[i][ii][k]=(float *)malloc(sizeof(****dp)*n[i]);
}
}
fprintf(stderr,"\nEnter %d x %d table of payoffs to the row player for state %d.\n",m[i],n[i],i+1);
for (k=0; k<m[i]; k++) {
for (j=0; j<n[i]; j++) {
if (scanf("%f",dp[i][0][k]+j)<1) badinp();
dp[i][0][k][j]/=scale;
}
}
for (ii=1; ii<=kk; ii++) {
fprintf(stderr,"\nEnter %d x %d transition table from state %d to state %d.\n",m[i],n[i],i+1,ii);
for (k=0; k<m[i]; k++) {
for (j=0; j<n[i]; j++) {
if (scanf("%f",dp[i][ii][k]+j)<1) badinp();
dp[i][ii][k][j]/=scale;
if (rec==3 && dp[i][ii][k][j]<0.) {
fprintf(stderr,"\nNegative transitions not allowed for total recursion.\n");
exit();
}
}
}
}
/*
* For total recursion make sure recursion adds up to one in every case.
*/
if (rec==3) {
for (k=0; k<m[i]; k++) {
for (j=0; j<n[i]; j++) {
v=0.;
for (ii=1; ii<=kk; ii++) v+=dp[i][ii][k][j];
if (v-1.<-1e-6 || v-1.>1e-6) {
fprintf(stderr,"\nRecursion of %f for state %d, row %d, column %d.\n",v,i+1,k+1,j+1);
exit();
}
}
}
}
}
}
/*
* Enter tables using sparse format.
*/
else {
for (i=0; i<nstate; i++) { /* allocate space for each state */
sp[i] = (struct entry ***)malloc(sizeof(**sp)*m[i]);
for (k=0; k<m[i]; k++) {
sp[i][k] = (struct entry **)malloc(sizeof(***sp)*n[i]);
for (j=0; j<n[i]; j++) sp[i][k][j] = NULL;
}
}
fprintf(stderr,"\nEnter payoff and transition table entries in the following formats:\n");
fprintf(stderr,"State Row Column 0 Payoff to row player\n");
fprintf(stderr,"State from Row Column State to Transition value\n");
fprintf(stderr,"0 (this signifies the end of table entries)\n");
while(1) {
if (scanf("%d",&i)<1) badinp();
if (i==0) break;
if (i<0 || i>nstate) {
fprintf(stderr,"\nBad state %d used for state from.\n",i);
exit();
}
if (scanf("%d%d%d%f",&k,&j,&ii,&v)<4) badinp();
if (k<1 || k>m[i-1]) {
fprintf(stderr,"\nBad row %d used for state %d.\n",k,i);
exit();
}
if (j<1 || j>n[i-1]) {
fprintf(stderr,"\nBad column %d used for state %d.\n",j,i);
exit();
}
if (ii<0 || ii>nstate || (rec==1 && ii>=i)) {
fprintf(stderr,"\nBad state %d used for state to.\n",ii);
exit();
}
if (v<0 && rec==3 && ii>0) {
fprintf(stderr,"\nNegative transition in total recursion mode.\n");
exit();
}
ent=sp[i-1][k-1][j-1];
while (ent) {
if (ent->estate==ii) {
fprintf(stderr,"\nDuplicate entry for state to %d, row %d, column %d, state from %d.\n",i,j,k,ii);
exit();
}
ent = ent->enext;
}
ent = (struct entry *)malloc(sizeof(*ent));
ent->enext=sp[i-1][k-1][j-1];
ent->eval=v/scale;
ent->estate=ii;
sp[i-1][k-1][j-1]=ent;
}
/*
* Make sure recursion adds up to one for total recursion.
*/
if (rec==3) {
for (i=0; i<nstate; i++) {
for (k=0; k<m[i]; k++) {
for (j=0; j<n[i]; j++) {
v=0.;
ent=sp[i][k][j];
while (ent) {
if (ent->estate!=0) v+=ent->eval;
ent=ent->enext;
}
if (v-1.<-1e6 || v-1.>1e-6) {
fprintf(stderr,"\nRecursion of %f for state %d, row %d, column %d.\n",v,i+1,k+1,j+1);
exit();
}
}
}
}
}
}
/*
* Allocate space for scratch matrix
*/
ii=nstate+1;
for (i=0; i<nstate; i++) {
if (m[i]>=ii) ii=m[i]+1;
}
a=(float **)malloc(sizeof(*a)*ii);
for (i=0; i<ii; i++) {
j=nstate+1;
for (k=0; k<nstate; k++) {
if (m[k]>=i && n[k]>=j) j=n[k]+1;
}
a[i]=(float *)malloc(sizeof(**a)*j);
}
/*
* Use zero for initial state values.
*/
old[0]=1.; /* value for payoff is special */
for (i=1; i<=nstate+1; i++) old[i]=0.;
/*
* Solve the game
*/
for (iter=1; ; iter++) { /* keep going until the answer is good enough */
/*
* Solve each state seperately.
*/
for (i=0; i<nstate; i++) {
if (rec==1) kk=i; /* recursion allowed? */
/*
* Build parameters for calling gamsol.
*/
if (inp==1) { /* dense input */
for (j=0; j<n[i]; j++) {
for (k=0; k<m[i]; k++) {
v=dp[i][0][k][j];
for (ii=1; ii<=kk; ii++) v+=old[ii]*dp[i][ii][k][j];
a[k][j]=v;
}
}
}
else { /* sparse input */
for (j=0; j<n[i]; j++) {
for (k=0; k<m[i]; k++) {
v=0.;
ent=sp[i][k][j];
while(ent) {
v+=ent->eval*old[ent->estate];
ent=ent->enext;
}
a[k][j]=v;
}
}
}
gamsol(m[i],n[i],a,row[i],col[i],val+i+1);
if (rec==1) old[i+1]=val[i+1]; /* if no recursion value can be used immediately */
}
/*
* Calculate state probabilities for total recursion.
* [ T ] [ ]
* [F'-I 1] X = [0] solved for X gives the state probabilities.
* [1 0] [1]
*/
if (rec==3) {
if (inp==1) { /* dense input */
for (i=0; i<nstate; i++) {
for (ii=0; ii<nstate; ii++) {
v=i==ii?-1.:0.;
for (j=0; j<n[i]; j++) {
for (k=0; k<m[i]; k++) {
v+=row[i][k]*col[i][j]*dp[i][ii+1][k][j];
}
}
a[ii][i]=v;
}
}
}
else { /* sparse input */
for (i=0; i<nstate; i++) {
for (ii=0; ii<nstate; ii++) a[ii][i]=i==ii?-1.:0.;
for (j=0; j<n[i]; j++) {
for (k=0; k<m[i]; k++) {
ent=sp[i][k][j];
while (ent) {
if (ent->estate!=0) a[ent->estate-1][i]+=row[i][k]*col[i][j]*ent->eval;
ent=ent->enext;
}
}
}
}
}
for (i=0; i<nstate; i++) {
b[i]=0;
a[i][nstate]=1.;
a[nstate][i]=1.;
}
b[nstate]=1.;
a[nstate][nstate]=0.;
i=lineqs(nstate+1,a,b,x);
if (i!=nstate+1) {
fprintf(stderr,"\nSingular matrix for probability calculation in iteration %d.\n",iter);
for (i=0; i<nstate; i++) x[i]=1./nstate;
}
val[nstate+1]=0.;
for (i=0; i<nstate; i++) val[nstate+1]+=x[i]*val[i+1];
for (i=1; i<=nstate; i++) val[i]-=old[nstate+1];
}
/*
* Check to see if we have done enough iterations.
*/
err=0.;
for (i=1; i<=nstate+1; i++) {
v=val[i]-old[i];
if (v<0.) v= -v;
if (v>err) err=v;
}
if (rec==1||iter>=miter||err<=merr) break;
/*
* Solve (I-F') X = Old - New, then set Old = Old - X. This will greatly speed
* up convergence if (I-F') is not singular. If it is then set Old = New.
*/
if (rec==2) { /* general recursion */
if (inp==1) { /* dense input */
for (i=0; i<nstate; i++) {
b[i]=old[i+1]-val[i+1];
for (ii=0; ii<nstate; ii++) {
v=i==ii?1.:0.;
for (j=0; j<n[i]; j++) {
for (k=0; k<m[i]; k++) {
v-=row[i][k]*col[i][j]*dp[i][ii+1][k][j];
}
}
a[i][ii]=v;
}
}
}
else { /* sparse input */
for (i=0; i<nstate; i++) {
b[i]=old[i+1]-val[i+1];
for (ii=0; ii<nstate; ii++) a[i][ii]=i==ii?1.:0.;
for (j=0; j<n[i]; j++) {
for (k=0; k<m[i]; k++) {
ent=sp[i][k][j];
while (ent) {
if (ent->estate!=0) a[i][ent->estate-1]-=row[i][k]*col[i][j]*ent->eval;
ent=ent->enext;
}
}
}
}
}
i=lineqs(nstate,a,b,x);
if (i==nstate) for (i=0; i<nstate; i++) old[i+1]=old[i+1]-x[i];
else {
fprintf(stderr,"\nSingular matrix in iteration %d.\n",iter);
for (i=1; i<=nstate; i++) old[i]=val[i];
}
}
/*
* ( [F -1])
* Solve (I - [P 0]) X = Old - New and then set Old = Old - X.
* This will greatly speed up convergence if the matrix is not singular.
* If it is set Old = New.
*/
else { /* total recursion */
if (inp==1) { /* dense input */
for (i=0; i<nstate; i++) {
for (ii=0; ii<nstate; ii++) {
v=i==ii?1.:0.;
for (j=0; j<n[i]; j++) {
for (k=0; k<m[i]; k++) {
v-=row[i][k]*col[i][j]*dp[i][ii+1][k][j];
}
}
a[i][ii]=v;
}
}
}
else { /* sparse input */
for (i=0; i<nstate; i++) {
for (ii=0; ii<nstate; ii++) a[i][ii]=i==ii?1.:0.;
for (j=0; j<n[i]; j++) {
for (k=0; k<m[i]; k++) {
ent=sp[i][k][j];
while (ent) {
if (ent->estate!=0) a[i][ent->estate-1]-=row[i][k]*col[i][j]*ent->eval;
ent=ent->enext;
}
}
}
}
}
for (i=0; i<=nstate; i++) b[i]=old[i+1]-val[i+1];
for (i=0; i<nstate; i++) {
a[i][nstate]=1.;
a[nstate][i]= -x[i]; /* copy state probabilities */
}
a[nstate][nstate]=1.;
i=lineqs(nstate+1,a,b,x);
if (i==nstate+1) for (i=0; i<=nstate; i++) old[i+1]=old[i+1]-x[i];
else {
fprintf(stderr,"\nSingular matrix in iteration %d.\n",iter);
for (i=1; i<=nstate; i++) old[i]=val[i]+old[nstate+1]-val[nstate+1];
old[nstate+1]=val[nstate+1];
}
}
}
/*
* Print out the results.
*/
if (rec!=1) {
printf("\nTermination after %d iterations.\n",iter);
printf("The largest value change in the last iteration was %e\n",err);
}
if (rec==3) {
printf("\nThe value per turn of the game is %.8f\n",val[nstate+1]);
printf("\nThe game state probabilities and values are:\n");
for (i=1; i<=nstate; i++) printf("State %4d %.8f %.8f\n",i,x[i-1],val[i]);
}
else {
printf("\nThe game state values are:\n");
for (i=1; i<=nstate; i++) printf("State %4d %.8f\n",i,val[i]);
}
for (i=0; i<nstate; i++) {
printf("\nThe active row strategies for state %d are:\n",i+1);
for (k=0; k<m[i]; k++) {
if (row[i][k]>=1.e-6) printf("Strategy %4d %.8f\n",k+1,row[i][k]);
}
printf("\nThe active column strategies for state %d are:\n",i+1);
for (j=0; j<n[i]; j++) {
if (col[i][j]>=1.e-6) printf("Strategy %4d %.8f\n",j+1,col[i][j]);
}
}
printf("\n");
}
/*
* Something about the input is messed up.
*/
void badinp()
{
fprintf(stderr,"\nBad input format or end of file reached.\n");
exit();
}
/*
* Subroutine to solve a zero sum game using the simplex method.
*
* Written by Bruno Wolff III
* June 13, 1989
*
* The game A is solved by minimizing the sum of X, with AX>=1, X>=0.
* For this to be solvable the value of the game must be positive. So
* that game A' is actually solved where A'=A-min(A)+1. The final answer
* is given by X=X/(sum of X). V=1/(sum of X)+min(A)-1. The dual solution
* is calculated at the same time and contains the other player's optimal
* strategy.
*
* Parameters:
* m Number of row strategies (input)
* n Number of column strategies (input)
* a Payoff matrix (to row player) with an extra row and column (input/scratch)
* r m optimal row strategies (output)
* c n optimal column strategies (output)
* v Value of the game to the row player (output)
*/
void gamsol(m,n,a,r,c,v)
int m;
int n;
float **a; /* a has an extra row and column to be used for work space */
float *r;
float *c;
float *v; /* pointer so as to allow a value to be returned */
{
float s; /* minimum value of the a array minus one */
float p; /* best value of the pivot criterion so far */
float pc; /* curent pivot criterion */
float cp; /* best value of the pivot criterion in the current column */
int ip,jp; /* location of best pivot */
int icp; /* location of the best pivot in the current column */
int i,j; /* scratch */
int *rin; /* active row strategies */
int *rout; /* inactive row strategies */
int *cin; /* active column strategies */
int *cout; /* inactive column strategies */
/*
* Check reasonableness of problem size.
*/
if (m<1) {
fprintf(stderr,"\ngamsol error - m must be at least 1.\n");
exit();
}
if (n<1) {
fprintf(stderr,"\ngamsol error - n must be at least 1.\n");
exit();
}
/*
* Make all payoffs positive. This assures the value of the game is
* positive.
*/
s=a[0][0]; /* start with s as some value in a */
for (i=0; i<m; i++) {
for (j=0; j<n; j++) {
if (a[i][j]<s) s=a[i][j];
}
}
s-=1.; /* this makes sure the value ends up positive */
for (i=0; i<m; i++) {
for (j=0; j<n; j++) {
a[i][j]-=s;
}
}
/*
* Finish setting up the simplex problem.
*/
rin=(int *)malloc(sizeof(*rin)*n);
rout=(int *)malloc(sizeof(*rout)*m);
cin=(int *)malloc(sizeof(*cin)*m);
cout=(int *)malloc(sizeof(*cout)*n);
for (i=0; i<m; i++) {
a[i][n]=1.;
rout[i]=i;
cin[i]= -1;
}
for (j=0; j<n; j++) {
a[m][j]= -1.;
cout[j]=j;
rin[j]= -1;
}
*v=1.; /* Value of the game times the sum of the strategies */
a[m][n]=0.;
/*
* Simplex loop
*/
while (1) { /* repeat until no suitable pivot is found */
/*
* Get pivot
*/
p= -1.; /* Mark that no pivot has been found yet */
for (j=0; j<n; j++) { /* Take the largest of the column values */
if (a[m][j]<0.) { /* Will a pivot in this column help? */
cp= -2.; /* Mark that no pivot has been found in this column */
for (i=0; i<m; i++) { /* Take the smallest value in each column */
if (a[i][j]>1.e-6) { /* Suitable pivots must be positive */
pc= -a[m][j]*a[i][n]/a[i][j];
if (pc<cp||cp== -2.) {
cp=pc;
icp=i;
}
}
}
if (cp>p) {
p=cp;
ip=icp;
jp=j;
}
}
}
/*
* If no pivot has been found the game is solved.
*/
if (p== -1.) break;
/*
* Perform the pivot operation.
*/
for (i=0; i<=m; i++) { /* Be sure to include extra row */
if (i!=ip) { /* Pivot row is unchanged except for the pivot */
for (j=0; j<=n; j++) { /* Be sure to include extra column */
if (j!=jp) { /* Negate pivot column (not the pivot) later */
a[i][j]=(a[i][j]*a[ip][jp]-a[ip][j]*a[i][jp])/(*v);
}
}
a[i][jp]= -a[i][jp]; /* Now it's safe to change pivot column */
}
}
p=a[ip][jp]; /* Don't update pivot until after we are done using it */
a[ip][jp]= *v; /* Switch pivot value with the value of the game */
*v=p;
i=rin[jp]; /* Keep track of active strategies */
rin[jp]=rout[ip];
rout[ip]=i;
j=cin[ip];
cin[ip]=cout[jp];
cout[jp]=j;
}
/*
* Copy over results.
*/
*v=(*v)/a[m][n]+s;
for (i=0; i<m; i++) r[i]=0; /* Set all row strategies to zero */
for (j=0; j<n; j++) { /* Then copy over the values of the active strategies */
if (rin[j]>=0) r[rin[j]]=a[m][j]/a[m][n];
}
for (j=0; j<n; j++) c[j]=0; /* Set all column strategies to zero */
for (i=0; i<m; i++) { /* Then copy over the values of the active strategies */
if (cin[i]>=0) c[cin[i]]=a[i][n]/a[m][n];
}
free(rin); /* don't let storage accumalate */
free(rout);
free(cin);
free(cout);
}
/*
* Function to solve AX=B by Guasian elimination.
*
* The rank of A is returned so that it is possible to tell when A is
* singular and X is not returned.
*
* Parameters:
* n The size of the problem (input)
* a n x n matrix (input/scratch)
* b n vector (input/scratch)
* x n vector (output)
*/
int lineqs(n,a,b,x)
int n;
float **a;
float *b;
float *x;
{
int *isw; /* Table used to keep track of row switches */
int *jsw; /* Table used to keep track of column switches */
int jp,kp; /* Pivot location */
float pivot; /* Pivot value */
int i,j,k; /* Scratch */
float e; /* Scratch */
/*
* Initialize row and column switch tables.
*/
isw=(int *)malloc(sizeof(*isw)*n);
jsw=(int *)malloc(sizeof(*jsw)*n);
for (i=0; i<n; i++) {
isw[i]=i;
jsw[i]=i;
}
/*
* Triangularize a.
*/
for (i=0; i<n; i++) {
/*
* Use maximal element as pivot.
*/
pivot=0.;
for (j=i; j<n; j++) {
for (k=i; k<n; k++) {
e=a[isw[j]][jsw[k]];
if (e<0.) e= -e;
if (e>pivot) {
pivot=e;
jp=j;
kp=k;
}
}
}
/*
* If the pivot is essentially zero then a is singular. If so we can't
* solve the problem and should just return the rank of a.
*/
if (pivot<1e-5) {
free(isw); /* don't let storage accumalate */
free(jsw);
return i;
}
/*
* Switch pivot row and column with the ith row and column.
*/
j=isw[i];
isw[i]=isw[jp];
isw[jp]=j;
k=jsw[i];
jsw[i]=jsw[kp];
jsw[kp]=k;
/*
* Eliminate the remainder of the column.
*/
pivot=a[isw[i]][jsw[i]];
for (j=i+1; j<n; j++) {
e= -a[isw[j]][jsw[i]]/pivot;
for (k=i+1; k<n; k++) {
a[isw[j]][jsw[k]]+=e*a[isw[i]][jsw[k]];
}
b[isw[j]]+=e*b[isw[i]];
}
}
/*
* Back substitute to obtain x.
*/
for (i=n-1; i>=0; i--) {
x[jsw[i]]=b[isw[i]];
for (j=n-1; j>i; j--) {
x[jsw[i]]-=a[isw[i]][jsw[j]]*x[jsw[j]];
}
x[jsw[i]]/=a[isw[i]][jsw[i]];
}
free(isw); /* don't let storage accumalate */
free(jsw);
return n; /* matrix is not singular, so return the full rank */
}
A secure version of this page is located at: https://wolff.to/bruno/strategy.html