Perform Gaussian elimination on the given matrix. This matrix has to be a
fat matrix (number of rows > number of columns).
435553
public void gaussianElimination(int[][] matrix) {
    assert(matrix != null && matrix.length > 0 && matrix[0].length > 0
           && matrix.length < matrix[0].length);
    int height = matrix.length;
    int width = matrix[0].length;
    for (int i = 0; i < height; i++) {
      boolean pivotFound = false;
      // scan the column for a nonzero pivot and swap it to the diagonal
      for (int j = i; j < height; j++) {
        if (matrix[i][j] != 0) {
          int[] tmp = matrix[i];
          matrix[i] = matrix[j];
          matrix[j] = tmp;
          pivotFound = true;
          break;
        }
      }
      if (!pivotFound) {
        continue;
      }
      int pivot = matrix[i][i];
      for (int j = i; j < width; j++) {
        matrix[i][j] = divide(matrix[i][j], pivot);
      }
      for (int j = i + 1; j < height; j++) {
        int lead = matrix[j][i];
        for (int k = i; k < width; k++) {
          matrix[j][k] = add(matrix[j][k], multiply(lead, matrix[i][k]));
        }
      }
    }
    for (int i = height - 1; i >=0; i--) {
      for (int j = 0; j < i; j++) {
        int lead = matrix[j][i];
        for (int k = i; k < width; k++) {
          matrix[j][k] = add(matrix[j][k], multiply(lead, matrix[i][k]));
        }
      }
    }
  }
64226
public static void doubleCenterMatrix(double[] dMatrix, int size) {
    double[] rowMean = new double[size];
    // row sum
    for(int i = 0, c = 0; i < size; i++) {
      for(int j = 0; j < i; j++) {
        double v = dMatrix[c++];
        rowMean[i] += v;
        rowMean[j] += v;
      }
      assert (dMatrix[c] == 0.);
      c++; // Diagonal entry. Must be zero!
    }
    // Normalize averages:
    double matrixMean = 0.;
    for(int i = 0; i < size; i++) {
      matrixMean += rowMean[i];
      rowMean[i] /= size;
    }
    matrixMean /= size * size;

    for(int o = 0, c = 0; o < size; o++) {
      // Including row mean!
      for(int p = 0; p <= o; p++) {
        dMatrix[c++] -= rowMean[o] + rowMean[p] - matrixMean;
      }
    }
  }
1240852
public static void byRow(Matrix m) {
        for (int i = 0; i < m.rows(); ++i) {
            double rowSum = 0;
            for (int j = 0; j < m.columns(); ++j)
                rowSum += m.get(i,j);

            // Skip row whose sume is zero.
            if (rowSum == 0) 
                continue;
            for (int j = 0; j < m.columns(); ++j)
                m.set(i,j, m.get(i,j) / rowSum);
        }
    }
67059
function dfFromMatrixHelper(matrix, startRow, colNames) {
  var nCol = colNames.size();
  var nRow = matrix.length - startRow;
  var columns = allocArray(nCol);
  var j;
  for (j = 0; j < nCol; j++) {
    columns[j] = allocArray(nRow);
  }
  for (var i = 0; i < nRow; i++) {
    var rowArray = matrix[i + startRow];
    if (rowArray.length !== nCol) {
      throw new Error('all row arrays must be of the same size');
    }
    for (j = 0; j < nCol; j++) {
      columns[j][i] = rowArray[j];
    }
  }
  return newDataFrame(columns, colNames);
}
1176828
public static void setMatrix(final double[][] m1, final int r0, final int r1, final int[] c, final double[][] m2) {
    assert r0 <= r1 : ERR_INVALID_RANGE;
    assert r1 <= m1.length : ERR_MATRIX_DIMENSIONS;
    for(int i = r0; i < r1; i++) {
      final double[] row1 = m1[i], row2 = m2[i - r0];
      for(int j = 0; j < c.length; j++) {
        row1[c[j]] = row2[j];
      }
    }
  }
1362057
public void decompose( DMatrixRMaj U , DMatrixRMaj S , DMatrixRMaj V ) {
		// this ensures the resulting rotation matrix will have a determinant of +1 and thus be a real rotation matrix
		if( CommonOps_DDRM.det(U) < 0 ) {
			CommonOps_DDRM.scale(-1,U);
			CommonOps_DDRM.scale(-1,S);
		}

		if( CommonOps_DDRM.det(V) < 0 ) {
			CommonOps_DDRM.scale(-1,V);
			CommonOps_DDRM.scale(-1,S);
		}

		// for possible solutions due to ambiguity in the sign of T and rotation
		extractTransform(U, V, S, solutions.get(0), true, true);
		extractTransform(U, V, S, solutions.get(1), true, false);
		extractTransform(U, V, S, solutions.get(2) , false,false);
		extractTransform(U, V, S, solutions.get(3), false, true);
	}
968104
public static void multAddColA(DMatrixSparseCSC A , int colA ,
                                   double alpha,
                                   DMatrixSparseCSC C, int mark,
                                   double x[] , int w[] ) {
        int idxA0 = A.col_idx[colA];
        int idxA1 = A.col_idx[colA+1];

        for (int j = idxA0; j < idxA1; j++) {
            int row = A.nz_rows[j];

            if( w[row] < mark ) {
                if( C.nz_length >= C.nz_rows.length ) {
                    C.growMaxLength(C.nz_length *2+1,true);
                }

                w[row] = mark;
                C.nz_rows[C.nz_length] = row;
                C.col_idx[mark] = ++C.nz_length;
                x[row] = A.nz_values[j]*alpha;
            } else {
                x[row] += A.nz_values[j]*alpha;
            }
        }
    }
340054
public boolean process( DMatrixSparseCSC A ) {
        init(A);

        TriangularSolver_DSCC.eliminationTree(A,true,parent,gwork);

        countNonZeroInR(parent);
        countNonZeroInV(parent);

        // if more columns than rows it's possible that Q*R != A. That's because a householder
        // would need to be created that's outside the  m by m Q matrix. In reality it has
        // a partial solution. Column pivot are needed.
        if( m < n ) {
            for (int row = 0; row <m; row++) {
                if( gwork.data[head+row] < 0 ) {
                    return false;
                }
            }
        }
        return true;
    }
533367
public static void mult( double a[] , int offsetA ,
                             DMatrixSparseCSC B ,
                             double c[] , int offsetC )
    {
        if( a.length-offsetA < B.numRows)
            throw new IllegalArgumentException("Length of 'a' isn't long enough");
        if( c.length-offsetC < B.numCols)
            throw new IllegalArgumentException("Length of 'c' isn't long enough");

        for (int k = 0; k < B.numCols; k++) {
            int idx0 = B.col_idx[k  ];
            int idx1 = B.col_idx[k+1];

            double sum = 0;
            for (int indexB = idx0; indexB < idx1; indexB++) {
                sum += a[offsetA+B.nz_rows[indexB]]*B.nz_values[indexB];
            }
            c[offsetC+k] = sum;
        }
    }
1123995
public void solve(DenseMatrix B) {
        if (!full) {
            throw new IllegalStateException("This is not a FULL singular value decomposition.");
        }

        if (B.nrows() != m) {
            throw new IllegalArgumentException("Dimensions do not agree.");
        }

        double[] b = new double[m];
        double[] x = new double[n];
        int p = B.ncols();
        for (int j = 0; j < p; j++) {
            for (int i = 0; i < m; i++) {
                b[i] = B.get(i, j);
            }

            solve(b, x);
            for (int i = 0; i < n; i++) {
                B.set(i, j, x[i]);
            }
        }
    }
1255453
def compute_affinity_matrix(adjacency_matrix, method='auto', **kwargs):
    """"""
    if method == 'auto':
        method = 'gaussian'
    return Affinity.init(method, **kwargs).affinity_matrix(adjacency_matrix)
542771
function(sbmA, sbmB) {
    // check if size compatible
    if (sbmA.y !== sbmB.x) {
      return
    }
    // create a new matrix
    return this.make((function (i, j) {
      return _.reduce(_.map(_.range(sbmA.y), (function(id) {
        return this.check(sbmA, i, id) && this.check(sbmB, id, j)
      }).bind(this)), function (sum, a) {
        return sum ^ a
      }) == 1
    }).bind(this), { x: sbmA.x, y: sbmB.y })
  }
784154
@Override
	public DoubleMatrix1D[] solve() throws Exception {

		final CholeskyFactorization HFact = new CholeskyFactorization(H, new Matrix1NornRescaler());
		boolean isHReducible = true;
		try{
			HFact.factorize();
		}catch(Exception e){
			isHReducible = false;
		}


        final DoubleMatrix1D v;// dim equals cols of A
        final DoubleMatrix1D w;// dim equals rank of A

		if (isHReducible) {
			// Solving KKT system via elimination
			DoubleMatrix1D HInvg = HFact.solve(g);
			if (A != null) {
				DoubleMatrix2D HInvAT = HFact.solve(AT);
				DoubleMatrix2D MenoSLower = ColtUtils.subdiagonalMultiply(A, HInvAT);
				DoubleMatrix1D AHInvg = ALG.mult(A, HInvg);
				
				CholeskyFactorization MSFact = new CholeskyFactorization(MenoSLower, new Matrix1NornRescaler());
				MSFact.factorize();
				if(h == null){
					w = MSFact.solve(ColtUtils.scalarMult(AHInvg, -1));
				}else{
					w = MSFact.solve(ColtUtils.add(h, AHInvg, -1));
				}
				
				v = HInvg.assign(ALG.mult(HInvAT, w), Functions.plus).assign(Mult.mult(-1));
			} else {
				w = null;
				v = HInvg.assign(Mult.mult(-1));
			}
		} else {
			// H is singular
			// Solving the full KKT system
			if(A!=null){
				DoubleMatrix1D[] fullSol =  this.solveAugmentedKKT();
				v = fullSol[0];
				w = fullSol[1];
			}else{
				//@TODO: try with rescaled H
				throw new Exception("KKT solution failed");
			}
		}

		// solution checking
		if (this.checkKKTSolutionAccuracy && !this.checkKKTSolutionAccuracy(v, w)) {
			log.error("KKT solution failed");
			throw new Exception("KKT solution failed");
		}

		DoubleMatrix1D[] ret = new DoubleMatrix1D[2];
		ret[0] = v;
		ret[1] = w;
		return ret;
	}
1145792
public static void multRow(Matrix A, int i, double c)
    {
        multRow(A, i, 0, A.cols(), c);
    }
368902
@Override
    public void solve(ZMatrixRMaj B, ZMatrixRMaj X) {
        if( X.numRows != numCols )
            throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numCols);
        else if( B.numRows != numRows || B.numCols != X.numCols )
            throw new IllegalArgumentException("Unexpected dimensions for B");

        U = decomposer.getR(U,true);
        final double gammas[] = decomposer.getGammas();
        final double dataQR[] = QR.data;

        final int BnumCols = B.numCols;

        // solve each column one by one
        for( int colB = 0; colB < BnumCols; colB++ ) {

            // make a copy of this column in the vector
            for( int i = 0; i < numRows; i++ ) {
                int indexB = (i*BnumCols + colB)*2;
                a[i*2]   = B.data[indexB];
                a[i*2+1] = B.data[indexB+1];
            }

            // Solve Qa=b
            // a = Q'b
            // a = Q_{n-1}...Q_2*Q_1*b
            //
            // Q_n*b = (I-gamma*u*u^H)*b = b - u*(gamma*U^H*b)
            for( int n = 0; n < numCols; n++ ) {
                int indexU = (n*numRows + n + 1)*2;

                double realUb = a[n*2];
                double imagUb = a[n*2+1];

                // U^H*b
                for( int i = n+1; i < numRows; i++ ) {
                    double realU = dataQR[indexU++];
                    double imagU = -dataQR[indexU++];

                    double realB = a[i*2];
                    double imagB = a[i*2+1];

                    realUb += realU*realB - imagU*imagB;
                    imagUb += realU*imagB + imagU*realB;
                }

                // gamma*U^T*b
                realUb *= gammas[n];
                imagUb *= gammas[n];

                a[n*2]   -= realUb;
                a[n*2+1] -= imagUb;

                indexU = (n*numRows + n + 1)*2;
                for( int i = n+1; i < numRows; i++) {
                    double realU = dataQR[indexU++];
                    double imagU = dataQR[indexU++];

                    a[i*2]   -= realU*realUb - imagU*imagUb;
                    a[i*2+1] -= realU*imagUb + imagU*realUb;
                }
            }

            // solve for Rx = b using the standard upper triangular solver
            TriangularSolver_ZDRM.solveU(U.data, a, numCols);

            // save the results

            for( int i = 0; i < numCols; i++ ) {
                int indexX = (i*X.numCols+colB)*2;

                X.data[indexX]   = a[i*2];
                X.data[indexX+1] = a[i*2+1];
            }
        }
    }
958542
@Override
    public Matrix[] decompose() {

        Matrix l = matrix.blankOfShape(matrix.rows(), matrix.rows());

        for (int j = 0; j < l.rows(); j++) {

            double d = 0.0;

            for (int k = 0; k < j; k++) {

                double s = 0.0;

                for (int i = 0; i < k; i++) {
                    s += l.get(k, i) * l.get(j, i);
                }

                s = (matrix.get(j, k) - s) / l.get(k, k);

                l.set(j, k, s);

                d = d + s * s;
            }

            d = matrix.get(j, j) - d;

            l.set(j, j, Math.sqrt(Math.max(d, 0.0)));

            for (int k = j + 1; k < l.rows(); k++) {
                l.set(j, k, 0.0);
            }
        }

        return new Matrix[] { l };
    }
395148
public static void emptyMatrix(byte[][] matrix, int maxX, int maxY) {
		for(int i = 0; i < maxX; i++) {
			for(int j = 0; j < maxY; j++) {
				matrix[i][j] = 0;
			}
		}
	}