OR the two matrices
19444
function(sbmA, sbmB) {
    // if not the same sizes
    if (sbmA.x != sbmB.x || sbmA.y != sbmB.y) {
      return
    }
    // create the matrix
    var mat = {}
    mat.x = sbmA.x
    mat.y = sbmA.y
    mat.data = []
    // populate
    for (var i = 0; i < mat.x; i++) {
      mat.data[i] = _.union(sbmA.data[i], sbmB.data[i]).sort()
    }
    // return matrix
    return mat
  }
1113302
function(a, b, r1, r2, k1, k2)
{
    // calculate mass matrix
    // If I wasn't lazy and wrote a proper matrix class, this wouldn't be so gross...
    var k11, k12, k21, k22;
    var m_sum = a.m_inv + b.m_inv;
    
    // start with I*m_sum
    k11 = m_sum; k12 = 0;
    k21 = 0;     k22 = m_sum;
    
    // add the influence from r1
    var a_i_inv = a.i_inv;
    var r1xsq =  r1.x * r1.x * a_i_inv;
    var r1ysq =  r1.y * r1.y * a_i_inv;
    var r1nxy = -r1.x * r1.y * a_i_inv;
    k11 += r1ysq; k12 += r1nxy;
    k21 += r1nxy; k22 += r1xsq;
    
    // add the influnce from r2
    var b_i_inv = b.i_inv;
    var r2xsq =  r2.x * r2.x * b_i_inv;
    var r2ysq =  r2.y * r2.y * b_i_inv;
    var r2nxy = -r2.x * r2.y * b_i_inv;
    k11 += r2ysq; k12 += r2nxy;
    k21 += r2nxy; k22 += r2xsq;
    
    // invert
    var determinant = k11*k22 - k12*k21;
    assertSoft(determinant !== 0, "Unsolvable constraint.");
    
    var det_inv = 1/determinant;

    k1.x =  k22*det_inv; k1.y = -k12*det_inv;
    k2.x = -k21*det_inv; k2.y =  k11*det_inv;
}
367879
public void multT(DoubleMatrix bMat, DoubleMatrix cMat, boolean transposeA, boolean transposeB) {
        if (!transposeA) {            
            throw new IllegalArgumentException("Multiplication where this matrix NOT transposed is not implemented.");
        }
        
        SparseRowDoubleMatrix.checkMultDimensions(this.viewTranspose(), transposeB ? bMat.viewTranspose() : bMat, cMat);

        int numRowsC = cMat.getNumRows();
        int numColsC = cMat.getNumColumns();
        
        if (bMat instanceof DenseDoubleMatrix) {
            DenseDoubleMatrix denseBMat = (DenseDoubleMatrix) bMat;
            for (int row = 0; row < numRowsC; row++) {
                for (int col = 0; col < numColsC; col++) {
                    if (transposeB) {
                        cMat.set(row, col, this.cols[row].dot(denseBMat.matrix[col]));
                    } else {
                        cMat.set(row, col, this.cols[row].dot(denseBMat.matrix, col));
                    }
                }
            }
        } else if (bMat instanceof SparseColDoubleMatrix && !transposeB) {
            SparseColDoubleMatrix sparseBMat = (SparseColDoubleMatrix) bMat;
            for (int row = 0; row < numRowsC; row++) {
                for (int col = 0; col < numColsC; col++) {
                    cMat.set(row, col, this.cols[row].dot(sparseBMat.cols[col]));
                }
            }
        } else if (bMat instanceof SparseRowDoubleMatrix && transposeB) {
            SparseRowDoubleMatrix sparseBMat = (SparseRowDoubleMatrix) bMat;
            for (int row = 0; row < numRowsC; row++) {
                for (int col = 0; col < numColsC; col++) {
                    cMat.set(row, col, this.cols[row].dot(sparseBMat.rows[col]));
                }
            }
        } else {
            throw new IllegalArgumentException("unhandled type: " + bMat.getClass().getCanonicalName());
        }
    }
211908
def concat(self, one, two):
        """"""
        if not len(one) == len(two) == 6:
            raise ValueError("bad sequ. length")
        self.a, self.b, self.c, self.d, self.e, self.f = TOOLS._concat_matrix(one, two)
        return self
798956
function _sizesMatch(matrix1, matrix2) {
  return matrix1.size.rows == matrix2.size.rows &&
    matrix1.size.columns == matrix2.size.columns;
}
1269769
public static DoubleMatrix solveSymmetric(DoubleMatrix A, DoubleMatrix B) {
		A.assertSquare();
		DoubleMatrix X = B.dup();
		int[] ipiv = new int[B.rows];
		SimpleBlas.sysv('U', A.dup(), ipiv, X);
		return X;
	}
225751
public static final DoubleMatrix2D add(DoubleMatrix2D A, DoubleMatrix2D B){
		if(A.rows()!=B.rows() || A.columns()!=B.columns()){
			throw new IllegalArgumentException("wrong matrices dimensions");
		}
		DoubleMatrix2D ret = DoubleFactory2D.dense.make(A.rows(), A.columns());
		for(int i=0; i<ret.rows(); i++){
			for(int j=0; j<ret.columns(); j++){
				ret.setQuick(i, j, A.getQuick(i, j) + B.getQuick(i, j));
			}
		}

		return ret;
	}
1367612
public static void multAddTransAB(DMatrixSparseCSC A , DMatrixRMaj B , DMatrixRMaj C )
    {
        if( A.numCols != C.numRows || B.numRows != C.numCols )
            throw new IllegalArgumentException("Inconsistent matrix shapes. "+stringShapes(A,B,C));

        ImplSparseSparseMult_DSCC.multAddTransAB(A,B,C);
    }
451433
public static void subtractEquals( DMatrix2 a , DMatrix2 b ) {
        a.a1 -= b.a1;
        a.a2 -= b.a2;
    }
633387
public static DenseMatrix block(Matrix a, Matrix b, Matrix c, Matrix d) {
        return Basic2DMatrix.block(a, b, c, d);
    }
1320237
function multiply_matrices(l, r) {
          var al = l[0];
          var cl = l[1];
          var el = l[2];
          var bl = l[3];
          var dl = l[4];
          var fl = l[5];

          var ar = r[0];
          var cr = r[1];
          var er = r[2];
          var br = r[3];
          var dr = r[4];
          var fr = r[5];

          var a = al * ar + cl * br;
          var c = al * cr + cl * dr;
          var e = al * er + cl * fr + el;
          var b = bl * ar + dl * br;
          var d = bl * cr + dl * dr;
          var f = bl * er + dl * fr + fl;

          return [a, c, e, b, d, f];
      }
46429
public static void addEquals( DMatrix2 a , DMatrix2 b ) {
        a.a1 += b.a1;
        a.a2 += b.a2;
    }
374400
public static void multTranA( DMatrixRMaj A , DMatrixRMaj B , DMatrixRMaj C , DMatrixRMaj output )
	{
		double t11 = A.data[0]*B.data[0] + A.data[3]*B.data[3] + A.data[6]*B.data[6];
		double t12 = A.data[0]*B.data[1] + A.data[3]*B.data[4] + A.data[6]*B.data[7];
		double t13 = A.data[0]*B.data[2] + A.data[3]*B.data[5] + A.data[6]*B.data[8];

		double t21 = A.data[1]*B.data[0] + A.data[4]*B.data[3] + A.data[7]*B.data[6];
		double t22 = A.data[1]*B.data[1] + A.data[4]*B.data[4] + A.data[7]*B.data[7];
		double t23 = A.data[1]*B.data[2] + A.data[4]*B.data[5] + A.data[7]*B.data[8];

		double t31 = A.data[2]*B.data[0] + A.data[5]*B.data[3] + A.data[8]*B.data[6];
		double t32 = A.data[2]*B.data[1] + A.data[5]*B.data[4] + A.data[8]*B.data[7];
		double t33 = A.data[2]*B.data[2] + A.data[5]*B.data[5] + A.data[8]*B.data[8];

		output.data[0] = t11*C.data[0] + t12*C.data[3] + t13*C.data[6];
		output.data[1] = t11*C.data[1] + t12*C.data[4] + t13*C.data[7];
		output.data[2] = t11*C.data[2] + t12*C.data[5] + t13*C.data[8];

		output.data[3] = t21*C.data[0] + t22*C.data[3] + t23*C.data[6];
		output.data[4] = t21*C.data[1] + t22*C.data[4] + t23*C.data[7];
		output.data[5] = t21*C.data[2] + t22*C.data[5] + t23*C.data[8];

		output.data[6] = t31*C.data[0] + t32*C.data[3] + t33*C.data[6];
		output.data[7] = t31*C.data[1] + t32*C.data[4] + t33*C.data[7];
		output.data[8] = t31*C.data[2] + t32*C.data[5] + t33*C.data[8];
	}
1212879
def linsolve(self, A, b):
        
        if self.sparselib == 'umfpack':
            return umfpack.linsolve(A, b)

        elif self.sparselib == 'klu':
            return klu.linsolve(A, b)
1275162
def _kp2(A, B):
    
    N = A.shape[0]
    if B.shape[0] != N:
        raise(ValueError)
    newshape1 = A.shape[1]*B.shape[1]
    return np.einsum('ijk,ilm->ijlkm', A, B).reshape(N, newshape1, -1)
244015
public static void checkMatrixMultiplication() {

        DoubleMatrix A = new DoubleMatrix(new double[][]{
                    {1.0, 2.0, 3.0},
                    {4.0, 5.0, 6.0},
                    {7.0, 8.0, 9.0}
                });
        DoubleMatrix E = new DoubleMatrix(new double[][]{
                    {0.0, 0.0, 1.0},
                    {0.0, 1.0, 0.0},
                    {1.0, 0.0, 0.0}
                });
        DoubleMatrix B = new DoubleMatrix(new double[][]{
                    {3.0, 2.0, 1.0},
                    {6.0, 5.0, 4.0},
                    {9.0, 8.0, 7.0}
                });

        check("checking matrix multiplication", A.mmul(E).equals(B));
    }
1086557
protected void eigenBB_2x2( int x1 ) {
        double b11 = diag[x1];
        double b12 = off[x1];
        double b22 = diag[x1+1];

        // normalize to reduce overflow
        double absA = Math.abs(b11);
        double absB = Math.abs(b12);
        double absC = Math.abs(b22);

        double scale = absA > absB ? absA : absB;
        if( absC > scale ) scale = absC;

        // see if it is a pathological case.  the diagonal must already be zero
        // and the eigenvalues are all zero.  so just return
        if( scale == 0 )
            return;

        b11 /= scale;
        b12 /= scale;
        b22 /= scale;

        eigenSmall.symm2x2_fast(b11*b11, b11*b12 , b12*b12+b22*b22);

        off[x1] = 0;
        diag[x1] = scale* Math.sqrt(eigenSmall.value0.real);
        double sgn = Math.signum(eigenSmall.value1.real);
        diag[x1+1] = sgn*scale* Math.sqrt(Math.abs(eigenSmall.value1.real));

    }