Reduced row echelon form: Difference between revisions

=={{header|D}}==

Be warned that this is definitely not the most efficient implementation, as it was coded from a thought process.

<lang d>

import std.stdio;

real[][]example = [

[1.0,2,-1,-4],

[2,3,-1,-11],

[-2,0,-3,22]

];

real getFirstNZ(real[]row) {

int i = getOrder(row);

if (i == row.length) return 0.0;

return row[i];

}

int getOrder(real[]row) {

foreach(i,ele;row) {

if (ele != 0) return i;

}

return row.length;

}

real[][]rref(real[][]input) {

real[][]ret;

// copy the input matrix

foreach(row;input) {

real[]currow = row.dup;

// fix leading digit to be positive to make things easier,

// by inverting the signs of all elements of the row if neccesary

if (getFirstNZ(currow) < 0) {

currow[] *= -1;

}

// normalize to the first element

currow[] /= getFirstNZ(currow);

ret ~= currow;

}

ret.sort.reverse;

// get the matrix into reduced row echelon form

for(int i = 0;i<ret.length;i++) {

int order = getOrder(ret[i]);

if (order == ret[i].length) {

// this row has no non-zero digits in it, so we're done,

// since any subsequent rows will also have all 0s, because of the sorting done earlier

break;

}

foreach(j,ref row;ret) {

// prevent us from erasing the current row

if (i == j) {

continue;

}

// see if we need to modify the current row

if (row[order] != 0.0) {

// cancel rows with awesome vector ops!

real[]tmp = ret[i].dup;

tmp[]*= -row[order];

row[]+=tmp[];

// ensure row is renormalized

int corder = getOrder(row);

if (corder < row.length && row[corder] != 1) {

row[]/=row[corder];

}

}

}

}

return ret;

}

void main() {

auto result = rref(example);

writefln(result);

}

</lang>