Inverting matrices

Ok, this is really annoying me.

I've attempted to do two different matrix inversion algorithms, but both of them destroy themselves on the first loop by dropping a zero on the leading diagonal and then trying to divide by it.

I know that this matrix can be inverted, so can anyone spot what I'm doing wrong? Please?

list matrix = [2.0,2.0,1.0,0.0,0.0, 2.0,2.0,2.0,1.0,0.0, 1.0,2.0,2.0,2.0,1.0, 0.0,1.0,2.0,2.0,2.0, 0.0,0.0,1.0,2.0,2.0];

integer size = 5;

float getMatrix(integer y, integer x) {
    return llList2Float(matrix, (y*size)+x);
}

setMatrix(integer y, integer x, float val) {
    integer addr = (y*size)+x;
    
    if (x == y) {
        if (val == 0.0) {
            llOwnerSay("Warning: Zero dropped on leading diagonal..");
        }
    }
    
    matrix = llListReplaceList(matrix,[val],addr,addr);
}


invertMatrix() {
    float swap;
    integer k;
    integer i;
    integer j;
    for (k=0; k<size; k++) {
        swap = getMatrix(k,k);
        setMatrix(k,k,1.0);
        for (i=0; i<size; i++) setMatrix(k,i,getMatrix(k,i) / swap);
        for (i=0; i<size; i++) if (i != k) {
            swap = getMatrix(i,k);
            setMatrix(i,k,0.0);
            for (j=0; j<size; j++) setMatrix(i,j,getMatrix(i,j) - (getMatrix(k,j) * swap));
        }
    }
}


colemanInvertMatrix() {
    integer k;
    integer i;
    integer j;
    
    for (k=0; k<size; k++) {
        llOwnerSay("k is " + (string)k);
        llSay(0, llDumpList2String(matrix,","));
               
        setMatrix(k,k,(-1/getMatrix(k,k)));
             
        for (i=0; i<size; i++) 
            if (i != k) setMatrix(i,k, getMatrix(i,k) * getMatrix(k,k) );
        
        for (i=0; i<size; i++) 
            if (i != k) 
                for (j=0; j<size; j++) 
                    if (j != k) {
                        setMatrix(i,j, getMatrix(i,j) + (getMatrix(i,k) * getMatrix(k,j)));
                    }
                    
        for (i=0; i<size; i++)
            if (i != k)
                setMatrix(k,i, getMatrix(k,i) * getMatrix(k,k));
    }
    
    for (i=0; i<size; i++) 
        for (j=0; j<size; j++) 
            setMatrix(i,j,(0-getMatrix(i,j)));
}
   

default
{
    
    touch_start(integer total_number)
    {
        llSay(0, llDumpList2String(matrix,","));
        invertMatrix();
        llSay(0, llDumpList2String(matrix,","));
    }
}

I haven't looked at the colemanInvertMatrix(), but I can see what's happening with the other one. The function is inverting the matrix using row operations, but it is ignoring the possibility of needing to swap rows in order to get a non zero entry (in the 2,2 spot in this case).

After using row one to kill off the rest of the first column, the 2,2 entry is zero. However the 3,2 entry (3rd row down, 2nd column) is non zero, and so in order to continue inverting you need to swap rows 2 and 3.

What seifert said, and also, I'm not sure how invertMatrix() is supposed to work anyway. You're reducing the matrix to the identity matrix, but that's not really going to help you here. The "standard" algorithm for finding a matrix's inverse is to tack on the identity matrix on the right side (in this case, a 5x5 identity). Then reduce the left 5x5 to the identity matrix, and what's left on the right will be the inverse. I'm not sure what the coleman inversion algorithm is, so I can't help you there.

As to actually reducing your matrix, it's faster if you only reduce coefficients to zero below the main diagonal first. Then start with the lower-rightmost leading 1, and use it to cancel all of the coefficients above it to zero. Continue on to the left with leading ones until you end up back at the top left again. I'm pretty sure that saves you multiplies, in the long run, which is really going to matter when implementing this in LSL.

I suggest looking up the Gauss-Jordan algorithm and implementing it to the letter. The algorithm includes catches for cases like Seifert's where a zero shows up and you have to swap rows up. Also, if you write code to hunt for a row to swap up, you'll also get the ability to see that there are no candidate rows with non-zero coefficients, which will mean that the matrix isn't invertible.

By the way, why're you doing linear algebra in SL? :)

Oh yeah, Lex is right, no augmented matrix stuff in there.

And linear algebra in R^5?