A symmetric matrix (M) is PD if it has a Cholesky decomposition, i.e.
    M = R.T dot R, where R is upper triangular with positive diagonal entries
559914
def check_pd(matrix):
    
    try:
        np.linalg.cholesky(matrix)
        return True
    except np.linalg.LinAlgError:
        return False
987540
public void factorize(boolean checkSymmetry) throws Exception {
		if (checkSymmetry && !Property.TWELVE.isSymmetric(Q)) {
			throw new Exception("Matrix is not symmetric");
		}

		if(this.rescaler != null){
			double[] cn_00_original = null;
			double[] cn_2_original = null;
			double[] cn_00_scaled = null;
			double[] cn_2_scaled = null;
			if(log.isDebugEnabled()){
				cn_00_original = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), Integer.MAX_VALUE);
				log.debug("cn_00_original Q before scaling: " + ArrayUtils.toString(cn_00_original));
				cn_2_original = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), 2);
				log.debug("cn_2_original Q before scaling : " + ArrayUtils.toString(cn_2_original));
			}
			//scaling the Q matrix, we have:
			//Q1 = U.Q.U[T] = U.L.L[T].U[T] = (U.L).(U.L)[T] 
			//and because U is diagonal it preserves the triangular form of U.L, so
			//Q1 = U.Q.U[T] = L1.L1[T] is the new Cholesky decomposition  
			DoubleMatrix1D Uv = rescaler.getMatrixScalingFactorsSymm(Q);
			if(log.isDebugEnabled()){
				boolean checkOK = rescaler.checkScaling(ColtUtils.fillSubdiagonalSymmetricMatrix(Q), Uv, Uv);
				if(!checkOK){
					log.warn("Scaling failed (checkScaling = false)");
				}
			}
			this.Q = ColtUtils.diagonalMatrixMult(Uv, Q, Uv);
			if(log.isDebugEnabled()){
				cn_00_scaled = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), Integer.MAX_VALUE);
				log.debug("cn_00_scaled Q after scaling : " + ArrayUtils.toString(cn_00_scaled));
				cn_2_scaled = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), 2);
				log.debug("cn_2_scaled Q after scaling  : " + ArrayUtils.toString(cn_2_scaled));
				
				if(cn_00_original[0] < cn_00_scaled[0] || cn_2_original[0] < cn_2_scaled[0]){
					//log.info("Q: " + ArrayUtils.toString(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()));
					log.warn("Problematic scaling");
					//throw new RuntimeException("Scaling failed");
				}
			}
		}

		switch (this.mode) {

		case DIAGONAL_PIVOLTING:
			pldltpt();
			break;

		case DIAGONAL_PIVOLTING_WITH_PARTIAL_PIVOTING:
			pldltptBK();
			break;

		default:
			throw new IllegalArgumentException("unknown mode " + mode);
		}
		
	}
71508
def inv_posdef(mat):
    """"""
    # See:
    # - http://www.seas.ucla.edu/~vandenbe/103/lectures/chol.pdf
    # - http://scicomp.stackexchange.com/questions/3188
    chol = np.linalg.cholesky(mat)
    ident = np.eye(mat.shape[0])
    res = solve_triangular(chol, ident, lower=True, overwrite_b=True)
    return np.transpose(res).dot(res)
219731
def convert_block_tridiag_to_banded(H_diag, H_upper_diag, lower=True):
    
    T, D, _ = H_diag.shape
    assert H_diag.ndim == 3 and H_diag.shape[2] == D
    assert H_upper_diag.shape == (T - 1, D, D)
    H_lower_diag = np.swapaxes(H_upper_diag, -2, -1)

    ab = np.zeros((2 * D, T * D))

    # Fill in blocks along the diagonal
    for d in range(D):
        # Get indices of (-d)-th diagonal of H_diag
        i = np.arange(d, D)
        j = np.arange(0, D - d)
        h = np.column_stack((H_diag[:, i, j], np.zeros((T, d))))
        ab[d] = h.ravel()

    # Fill in lower left corner of blocks below the diagonal
    for d in range(0, D):
        # Get indices of (-d)-th diagonal of H_diag
        i = np.arange(d, D)
        j = np.arange(0, D - d)
        h = np.column_stack((H_lower_diag[:, i, j], np.zeros((T - 1, d))))
        ab[D + d, :D * (T - 1)] = h.ravel()

    # Fill in upper corner of blocks below the diagonal
    for d in range(1, D):
        # Get indices of (+d)-th diagonal of H_lower_diag
        i = np.arange(0, D - d)
        j = np.arange(d, D)
        h = np.column_stack((np.zeros((T - 1, d)), H_lower_diag[:, i, j]))
        ab[D - d, :D * (T - 1)] += h.ravel()

    return ab if lower else transpose_lower_banded_matrix(ab)
723425
def _inverse_lower_triangular(M):
    
    if M.get_shape().ndims != 3:  # pragma: no cover
        raise ValueError("Number of dimensions for input is required to be 3.")
    D, N = tf.shape(M)[0], tf.shape(M)[1]
    I_DNN = tf.eye(N, dtype=M.dtype)[None, :, :] * tf.ones((D, 1, 1), dtype=M.dtype)
    return tf.matrix_triangular_solve(M, I_DNN)
1335003
public static INDArray symmetricGeneralizedEigenvalues(INDArray A, INDArray B) {
        Preconditions.checkArgument(A.isMatrix() && A.isSquare(), "Argument A must be a square matrix: has shape %s", A.shape());
        Preconditions.checkArgument(B.isMatrix() && B.isSquare(), "Argument B must be a square matrix: has shape %s", B.shape());
        INDArray W = Nd4j.create(A.rows());

        A = InvertMatrix.invert(B, false).mmuli(A);
        Nd4j.getBlasWrapper().syev('V', 'L', A, W);
        return W;
    }
1259653
private Intrinsic solveForCalibration(DMatrixRMaj w) {
		Intrinsic calib = new Intrinsic();

//		CholeskyDecomposition_F64<DMatrixRMaj> chol = DecompositionFactory_DDRM.chol(false);
//
//		chol.decompose(w.copy());
//		DMatrixRMaj R = chol.getT(w);
//		R.print();

		if( zeroSkew ) {
			calib.skew = 0;
			calib.fy = Math.sqrt(w.get(1,1));

			if( knownAspect ) {
				calib.fx = calib.fy/aspectRatio;
			} else {
				calib.fx = Math.sqrt(w.get(0,0));
			}
		} else if( knownAspect ) {
			calib.fy = Math.sqrt(w.get(1,1));
			calib.fx = calib.fy/aspectRatio;
			calib.skew = w.get(0,1)/calib.fy;
		} else {
			calib.fy = Math.sqrt(w.get(1,1));
			calib.skew = w.get(0,1)/calib.fy;
			calib.fx = Math.sqrt(w.get(0,0) - calib.skew*calib.skew);
		}
		return calib;
	}
1067945
def cublasDsymv(handle, uplo, n, alpha, A, lda, x, incx, beta, y, incy):
    

    status = _libcublas.cublasDsymv_v2(handle, 
                                       _CUBLAS_FILL_MODE[uplo], n, 
                                       ctypes.byref(ctypes.c_double(alpha)), 
                                       int(A), lda, int(x), incx, 
                                       ctypes.byref(ctypes.c_double(beta)), 
                                       int(y), incy)
    cublasCheckStatus(status)
476306
public Matrix4x3d lerp(Matrix4x3dc other, double t, Matrix4x3d dest) {
        dest.m00 = m00 + (other.m00() - m00) * t;
        dest.m01 = m01 + (other.m01() - m01) * t;
        dest.m02 = m02 + (other.m02() - m02) * t;
        dest.m10 = m10 + (other.m10() - m10) * t;
        dest.m11 = m11 + (other.m11() - m11) * t;
        dest.m12 = m12 + (other.m12() - m12) * t;
        dest.m20 = m20 + (other.m20() - m20) * t;
        dest.m21 = m21 + (other.m21() - m21) * t;
        dest.m22 = m22 + (other.m22() - m22) * t;
        dest.m30 = m30 + (other.m30() - m30) * t;
        dest.m31 = m31 + (other.m31() - m31) * t;
        dest.m32 = m32 + (other.m32() - m32) * t;
        dest.properties = properties & other.properties();
        return dest;
    }
403717
function d3_geom_hullCCW(i1, i2, i3, v) {
  var t, a, b, c, d, e, f;
  t = v[i1]; a = t[0]; b = t[1];
  t = v[i2]; c = t[0]; d = t[1];
  t = v[i3]; e = t[0]; f = t[1];
  return ((f-b)*(c-a) - (d-b)*(e-a)) > 0;
}
30541
def _embedding_classical_mds(matrix, dimensions=3, additive_correct=False):
    
    if additive_correct:
        dbc = double_centre(_additive_correct(matrix))
    else:
        dbc = double_centre(matrix)
    decomp = eigen(dbc)
    lambda_ = np.diag(np.sqrt(np.abs(decomp.vals[:dimensions])))
    evecs = decomp.vecs[:, :dimensions]
    coords = evecs.dot(lambda_)
    return coords
1254670
public static FloatMatrix symmetricGeneralizedEigenvalues(FloatMatrix A, FloatMatrix B) {
        A.assertSquare();
        B.assertSquare();
        FloatMatrix W = new FloatMatrix(A.rows);
        SimpleBlas.sygvd(1, 'N', 'U', A.dup(), B.dup(), W);
        return W;
    }
1337444
def issymmetric(am):
    
    """
    if am.shape[0] != am.shape[1]:
        logger.warning('Matrix is not square, symmetrical is irrelevant')
        return False
    if am.format != 'coo':
        am = am.tocoo(copy=False)
    if istril(am) or istriu(am):
        return False
    # Compare am with its transpose, element wise
    sym = ((am != am.T).size) == 0
    return sym
1241915
public Matrix4d lerp(Matrix4dc other, double t, Matrix4d dest) {
        dest.m00 = m00 + (other.m00() - m00) * t;
        dest.m01 = m01 + (other.m01() - m01) * t;
        dest.m02 = m02 + (other.m02() - m02) * t;
        dest.m03 = m03 + (other.m03() - m03) * t;
        dest.m10 = m10 + (other.m10() - m10) * t;
        dest.m11 = m11 + (other.m11() - m11) * t;
        dest.m12 = m12 + (other.m12() - m12) * t;
        dest.m13 = m13 + (other.m13() - m13) * t;
        dest.m20 = m20 + (other.m20() - m20) * t;
        dest.m21 = m21 + (other.m21() - m21) * t;
        dest.m22 = m22 + (other.m22() - m22) * t;
        dest.m23 = m23 + (other.m23() - m23) * t;
        dest.m30 = m30 + (other.m30() - m30) * t;
        dest.m31 = m31 + (other.m31() - m31) * t;
        dest.m32 = m32 + (other.m32() - m32) * t;
        dest.m33 = m33 + (other.m33() - m33) * t;
        dest.properties = properties & other.properties();
        return dest;
    }
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);
	}
341621
def cublasSdgmm(handle, mode, m, n, A, lda, x, incx, C, ldc):
    

    status = _libcublas.cublasSdgmm(handle,
                                    _CUBLAS_SIDE[mode],
                                    m, n, 
                                    int(A), lda, 
                                    int(x), incx,
                                    int(C), ldc)
    cublasCheckStatus(status)
822504
public static FloatMatrix symmetricEigenvalues(FloatMatrix A) {
        A.assertSquare();
        FloatMatrix eigenvalues = new FloatMatrix(A.rows);
        int isuppz[] = new int[2 * A.rows];
        SimpleBlas.syevr('N', 'A', 'U', A.dup(), 0, 0, 0, 0, 0, eigenvalues, dummyFloat, isuppz);
        return eigenvalues;
    }