I reedited this post, I obviously did not implement ludcmp well, apologies for that. Here's a function that does work, at least for the test case I tried, in case anyone is interested

int GenMat::ludcmp_nr(int indx[]) {
if (GetDim()!=2 || GetExt(0)!=GetExt(1)) {
cout<<"Matrix cannot be diagonalized in ludcmp\n";
exit(0);
}
unsigned int imax,n=GetExt(0);
float big,dum,vv[n];
float sum;
int inv=1; //No row interchanges yet.
for (unsigned int i=0;i<n;i++) {
big=0.0;float temp;
for (unsigned int j=0;j<n;j++)
if ((temp=fabs(m(Index(i,j)))) > big) big=temp;
if (big == 0.0) exit(0);
vv[i]=1.0/big;
}
for (unsigned int j=0;j<n;j++) {
int jj=Index(j,j);
for (unsigned int i=0;i<j;i++) {
sum=m(Index(i,j));
for (unsigned int k=0;k<i;k++) sum -= m(Index(i,k))*m(Index(k,j));
Setm(Index(i,j),sum);
}
big=0.0;
for (unsigned int i=j;i<n;i++) {
int ij=Index(i,j);
float sum=m(ij);
for (unsigned int k=0;k<j;k++) sum -=m(Index(i,k))*m(Index(k,j));
Setm(ij,sum);
if ((dum=vv[i]*fabs(sum)) >= big) {
//Is the figure of merit for the pivot better than the best so far?
big=dum;
imax=i;
}
}
if (j != imax) { // Do we need to interchange rows?
for (unsigned int k=0;k<n;k++) { // Yes, do so...
int imk=Index(imax,k); int jk=Index(j,k);
dum=m(imk);
Setm(imk,m(jk));
Setm(jk,dum);
}
inv = -(inv); // ...and change the parity of inv
vv[imax]=vv[j];
}
indx[j]=imax;
if (m(jj) != 0.0 && j!=n-1) {
dum=1/m(jj);
for (unsigned int i=j+1;i<n;i++) {
int ij=Index(i,j);
Mulm(ij,dum);
}
}
}
return inv;
}

Cheers,
Roeland.