// -----------------------------------------------------------------------------------
function gauss_jordan_solve(matrix) {   // Solve a system of simultaneous linear equations
                    // We'll just do Gauss-Jordan elimination

var nrows = matrix.length;
var ncols = matrix[0].length;   // ncols should be nrows+1, assuming it's a square matrix
 
var temprow = [];  // A temp for swapping two rows
var maxd;     // The max value in the diag-th column, on or below the pivotrow-th row
var maxdrow;  // The row that maxd is in
var tolerance = .0001;  // Any pivot less than this is considered to be zero
var pivots = [];
var pivotrow = 0;
var pivot = 0.3

// alert("nrows, ncols:  " + nrows + "   " + ncols + "\n");
// dump2D(matrix); 
for (var diag=0; diag<ncols-1; diag++) {   // For each column ...
        var i;

        pivots[diag] = -1;
        
        // Find maxd, the max absolute value in the diag-th column, on or below the pivotrow-th row
        // alert("Looking for maxd, for diag = " + diag + "\n");
        maxd = Math.abs(matrix[pivotrow][diag]);
        maxdrow = pivotrow;
        for (i=pivotrow+1; i<nrows; i++) {
            if (Math.abs(matrix[i][diag]) > maxd ){
                maxd = Math.abs(matrix[i][diag]);
                maxdrow = i;
            }
        }  
        // alert("diag = " + diag + "\n");
        //alert("maxd = " + maxd + "\n");
        //alert("maxdrow = " + maxdrow + "\n");
        
        // Now maxd and maxdrow have been set.
        // If pivotrow and maxdrow are not equal, swap the maxdrow-th row with the pivotrow-th row:
        if (pivotrow < maxdrow) {  // Swap rows
            // alert("Swapping rows " + diag + " and " + maxdrow + "\n");
            for (i =0; i<ncols; i++) { temprow[i] = matrix[pivotrow][i]; }          // Stash row pivotrow into temprow[]
            for (i =0; i<ncols; i++) { matrix[pivotrow][i] = matrix[maxdrow][i]; }  // Copy row maxdrow to row diag
            for (i =0; i<ncols; i++) { matrix[maxdrow][i] = temprow[i]; }       // Copy temprow[] to row maxdrow
            // alert("Swapped rows\n");
        }
        

        pivot = matrix[pivotrow][diag];
        // alert("Ready to pivot.  pivot = " + pivot + " tolerance = " + tolerance + "\n");
        // Now the pivotrow-th row has the pivot
        
        
        if (Math.abs(pivot) > tolerance) {  // pivot
        
            // alert("diag, pivotrow: " + diag + " " + pivotrow + "\n");
            pivots[diag] = pivotrow;

            for (i=0; i<ncols; i++) {matrix[pivotrow][i] /= pivot;}  // Put a 1 on the diagonal at matrix[pivotrow][diag]
            
            // Add the correct multiple of matrix[pivotrow] to each other row to zero out the diag-th column
            for (var irow=0; irow<nrows; irow++) {
                if (irow == pivotrow) {continue;}
                var fac = matrix[irow][diag];
                for (i=0; i<ncols; i++) {matrix[irow][i] = matrix[irow][i] - matrix[pivotrow][i]*fac;}    
            }
            pivotrow++;
        }
        else  {      // if (pivot < tolerance) just go to the next column of the matrix
            // alert("pivot is zero.  diag, pivotrow: " + diag + " " + pivotrow + "\n");
        }
}

// The matrix is now in reduced row echelon form: zeroes above and below each pivot,
//   and each pivot equal to 1.

// alert("nrows, ncols, pivotrow: " + nrows + " " + ncols + " " + pivotrow + "\n");
// alert("pivots = " + pivots + "\n");
// dump2D(matrix);

//  Use the pivots[] information to move the strengths to the right places:

for (i=nrows-1; i>=0; i--) {
    if (pivots[i] != -1) { 
        matrix[i][ncols-1] = matrix[pivots[i]][ncols-1]; 
    }
    else {
        matrix[i][ncols-1] = 0;
        // alert("setting i to zero:" + i + "\n");
    }
}  



return matrix;  
   
}     // End of solve()