A -> [L,U] A=LU L is lower triangular matrix U is upper triangular matrix
480990
function lu(A) {
    var size = A.length;
    //var L=jStat.diagonal(jStat.ones(1,size)[0]);
    var L = jStat.identity(size);
    var R = jStat.zeros(A.length, A[0].length);
    var parts;
    jStat.arange(size).forEach(function(t) {
      R[0][t] = A[0][t];
    });
    jStat.arange(1, size).forEach(function(l) {
      jStat.arange(l).forEach(function(i) {
        parts = jStat.arange(i).map(function(jj) {
          return L[l][jj] * R[jj][i];
        });
        L[l][i] = (A[l][i] - jStat.sum(parts)) / R[i][i];
      });
      jStat.arange(l, size).forEach(function(j) {
        parts = jStat.arange(l).map(function(jj) {
          return L[l][jj] * R[jj][j];
        });
        R[l][j] = A[i][j] - jStat.sum(parts);
      });
    });
    return [L, R];
  }
726928
public static boolean cholL( DMatrix6x6 A ) {

        A.a11 = Math.sqrt(A.a11);
        A.a12 = 0;
        A.a13 = 0;
        A.a14 = 0;
        A.a15 = 0;
        A.a16 = 0;
        A.a21 = (A.a21)/A.a11;
        A.a22 = Math.sqrt(A.a22-A.a21*A.a21);
        A.a23 = 0;
        A.a24 = 0;
        A.a25 = 0;
        A.a26 = 0;
        A.a31 = (A.a31)/A.a11;
        A.a32 = (A.a32-A.a31*A.a21)/A.a22;
        A.a33 = Math.sqrt(A.a33-A.a31*A.a31-A.a32*A.a32);
        A.a34 = 0;
        A.a35 = 0;
        A.a36 = 0;
        A.a41 = (A.a41)/A.a11;
        A.a42 = (A.a42-A.a41*A.a21)/A.a22;
        A.a43 = (A.a43-A.a41*A.a31-A.a42*A.a32)/A.a33;
        A.a44 = Math.sqrt(A.a44-A.a41*A.a41-A.a42*A.a42-A.a43*A.a43);
        A.a45 = 0;
        A.a46 = 0;
        A.a51 = (A.a51)/A.a11;
        A.a52 = (A.a52-A.a51*A.a21)/A.a22;
        A.a53 = (A.a53-A.a51*A.a31-A.a52*A.a32)/A.a33;
        A.a54 = (A.a54-A.a51*A.a41-A.a52*A.a42-A.a53*A.a43)/A.a44;
        A.a55 = Math.sqrt(A.a55-A.a51*A.a51-A.a52*A.a52-A.a53*A.a53-A.a54*A.a54);
        A.a56 = 0;
        A.a61 = (A.a61)/A.a11;
        A.a62 = (A.a62-A.a61*A.a21)/A.a22;
        A.a63 = (A.a63-A.a61*A.a31-A.a62*A.a32)/A.a33;
        A.a64 = (A.a64-A.a61*A.a41-A.a62*A.a42-A.a63*A.a43)/A.a44;
        A.a65 = (A.a65-A.a61*A.a51-A.a62*A.a52-A.a63*A.a53-A.a64*A.a54)/A.a55;
        A.a66 = Math.sqrt(A.a66-A.a61*A.a61-A.a62*A.a62-A.a63*A.a63-A.a64*A.a64-A.a65*A.a65);
        return !UtilEjml.isUncountable(A.a66);
    }
1165072
def restricted_cover(l, succsOf):
    
    fzl = frozenset(l)
    lut = dict()
    for i in l:
        lut[i] = fzl.intersection(succsOf(i))
    return lambda x: lut[x]
898384
def _cho_solve_AATI(A, rho, b, c, lwr, check_finite=True):
    """"""

    N, M = A.shape
    if N >= M:
        x = (b - _cho_solve((c, lwr), b.dot(A).T,
            check_finite=check_finite).T.dot(A.T)) / rho
    else:
        x = _cho_solve((c, lwr), b.T, check_finite=check_finite).T
    return x
1035019
private static int[] elmhes(DenseMatrix A) {
        int n = A.nrows();
        int[] perm = new int[n];

        for (int m = 1; m < n - 1; m++) {
            double x = 0.0;
            int i = m;
            for (int j = m; j < n; j++) {
                if (Math.abs(A.get(j, m - 1)) > Math.abs(x)) {
                    x = A.get(j, m - 1);
                    i = j;
                }
            }
            perm[m] = i;
            if (i != m) {
                for (int j = m - 1; j < n; j++) {
                    double swap = A.get(i, j);
                    A.set(i, j, A.get(m, j));
                    A.set(m, j, swap);
                }
                for (int j = 0; j < n; j++) {
                    double swap = A.get(j, i);
                    A.set(j, i, A.get(j, m));
                    A.set(j, m, swap);
                }
            }
            if (x != 0.0) {
                for (i = m + 1; i < n; i++) {
                    double y = A.get(i, m - 1);
                    if (y != 0.0) {
                        y /= x;
                        A.set(i, m - 1, y);
                        for (int j = m; j < n; j++) {
                            A.sub(i, j, y * A.get(m, j));
                        }
                        for (int j = 0; j < n; j++) {
                            A.add(j, m, y * A.get(j, i));
                        }
                    }
                }
            }
        }

        return perm;
    }
1373115
def cublasDtrmv(handle, uplo, trans, diag, n, A, lda, x, inx):
    

    status = _libcublas.cublasDtrmv_v2(handle, 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       n, int(A), lda, int(x), inx)
    cublasCheckStatus(status)
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)
1230368
public static void jacobian_Control3( DMatrixRMaj L_full ,
										  double beta[] ,  DMatrixRMaj A)
	{
		int indexA = 0;

		double b0 = beta[0]; double b1 = beta[1]; double b2 = beta[2];

		final double ld[] = L_full.data;

		for( int i = 0; i < 3; i++ ) {
			int li = L_full.numCols*i;
			A.data[indexA++] = 2*ld[li+0]*b0 +   ld[li+1]*b1 +   ld[li+2]*b2;
			A.data[indexA++] =   ld[li+1]*b0 + 2*ld[li+3]*b1 +   ld[li+4]*b2;
			A.data[indexA++] =   ld[li+2]*b0 +   ld[li+4]*b1 + 2*ld[li+5]*b2;
		}
	}
1186035
public static void multAddOuter( double alpha , DMatrix6x6 A , double beta , DMatrix6 u , DMatrix6 v , DMatrix6x6 C ) {
        C.a11 = alpha*A.a11 + beta*u.a1*v.a1;
        C.a12 = alpha*A.a12 + beta*u.a1*v.a2;
        C.a13 = alpha*A.a13 + beta*u.a1*v.a3;
        C.a14 = alpha*A.a14 + beta*u.a1*v.a4;
        C.a15 = alpha*A.a15 + beta*u.a1*v.a5;
        C.a16 = alpha*A.a16 + beta*u.a1*v.a6;
        C.a21 = alpha*A.a21 + beta*u.a2*v.a1;
        C.a22 = alpha*A.a22 + beta*u.a2*v.a2;
        C.a23 = alpha*A.a23 + beta*u.a2*v.a3;
        C.a24 = alpha*A.a24 + beta*u.a2*v.a4;
        C.a25 = alpha*A.a25 + beta*u.a2*v.a5;
        C.a26 = alpha*A.a26 + beta*u.a2*v.a6;
        C.a31 = alpha*A.a31 + beta*u.a3*v.a1;
        C.a32 = alpha*A.a32 + beta*u.a3*v.a2;
        C.a33 = alpha*A.a33 + beta*u.a3*v.a3;
        C.a34 = alpha*A.a34 + beta*u.a3*v.a4;
        C.a35 = alpha*A.a35 + beta*u.a3*v.a5;
        C.a36 = alpha*A.a36 + beta*u.a3*v.a6;
        C.a41 = alpha*A.a41 + beta*u.a4*v.a1;
        C.a42 = alpha*A.a42 + beta*u.a4*v.a2;
        C.a43 = alpha*A.a43 + beta*u.a4*v.a3;
        C.a44 = alpha*A.a44 + beta*u.a4*v.a4;
        C.a45 = alpha*A.a45 + beta*u.a4*v.a5;
        C.a46 = alpha*A.a46 + beta*u.a4*v.a6;
        C.a51 = alpha*A.a51 + beta*u.a5*v.a1;
        C.a52 = alpha*A.a52 + beta*u.a5*v.a2;
        C.a53 = alpha*A.a53 + beta*u.a5*v.a3;
        C.a54 = alpha*A.a54 + beta*u.a5*v.a4;
        C.a55 = alpha*A.a55 + beta*u.a5*v.a5;
        C.a56 = alpha*A.a56 + beta*u.a5*v.a6;
        C.a61 = alpha*A.a61 + beta*u.a6*v.a1;
        C.a62 = alpha*A.a62 + beta*u.a6*v.a2;
        C.a63 = alpha*A.a63 + beta*u.a6*v.a3;
        C.a64 = alpha*A.a64 + beta*u.a6*v.a4;
        C.a65 = alpha*A.a65 + beta*u.a6*v.a5;
        C.a66 = alpha*A.a66 + beta*u.a6*v.a6;
    }
278387
def matrix2lha(M):
    """"""
    l = []
    ind = np.indices(M.shape).reshape(M.ndim, M.size).T
    for i in ind:
        l.append([j+1 for j in i] + [M[tuple(i)]])
    return l
892904
def split_linear_constraints(A, l, u):
    
    ieq = []
    igt = []
    ilt = []
    ibx = []
    for i in range(len(l)):
        if abs(u[i] - l[i]) <= EPS:
            ieq.append(i)
        elif (u[i] > 1e10) and (l[i] > -1e10):
            igt.append(i)
        elif (l[i] <= -1e10) and (u[i] < 1e10):
            ilt.append(i)
        elif (abs(u[i] - l[i]) > EPS) and (u[i] < 1e10) and (l[i] > -1e10):
            ibx.append(i)
        else:
            raise ValueError

    Ae = A[ieq, :]
    Ai = sparse([A[ilt, :], -A[igt, :], A[ibx, :], -A[ibx, :]])
    be = u[ieq, :]
    bi = matrix([u[ilt], -l[igt], u[ibx], -l[ibx]])

    return Ae, be, Ai, bi
420224
def cublasDsyr(handle, uplo, n, alpha, x, incx, A, lda):
    

    status = _libcublas.cublasDsyr_v2(handle,
                                      _CUBLAS_FILL_MODE[uplo], n, 
                                      ctypes.byref(ctypes.c_double(alpha)), 
                                      int(x), incx, int(A), lda)
    cublasCheckStatus(status)
1321172
private void fillKsubMatrix( Coordinate[] controlPoints, GeneralMatrix L ) {
        double alfa = 0;
        int controlPointsNum = controlPoints.length;
        for( int i = 0; i < controlPointsNum; i++ ) {
            for( int j = i + 1; j < controlPointsNum; j++ ) {
                double u = calculateFunctionU(controlPoints[i], controlPoints[j]);
                L.setElement(i, j, u);
                L.setElement(j, i, u);
                alfa = alfa + (u * 2); // same for upper and lower part
            }
        }

        alfa = alfa / (controlPointsNum * controlPointsNum);
    }
622021
def dpotrs(A, B, lower=1):
    
    A = force_F_ordered(A)
    return lapack.dpotrs(A, B, lower=lower)
305442
public static FloatMatrix cholesky(FloatMatrix A) {
        FloatMatrix result = A.dup();
        int info = NativeBlas.spotrf('U', A.rows, result.data, 0, A.rows);
        if (info < 0) {
            throw new LapackArgumentException("DPOTRF", -info);
        } else if (info > 0) {
            throw new LapackPositivityException("DPOTRF", "Minor " + info + " was negative. Matrix must be positive definite.");
        }
        clearLower(result);
        return result;
    }
827776
def aorb(a, b):
    """"""
    return matrix(np.logical_or(a, b).astype('float'), a.size)
545209
private static void eltran(DenseMatrix A, DenseMatrix V, int[] perm) {
        int n = A.nrows();
        for (int mp = n - 2; mp > 0; mp--) {
            for (int k = mp + 1; k < n; k++) {
                V.set(k, mp, A.get(k, mp - 1));
            }
            int i = perm[mp];
            if (i != mp) {
                for (int j = mp; j < n; j++) {
                    V.set(mp, j, V.get(i, j));
                    V.set(i, j, 0.0);
                }
                V.set(i, mp, 1.0);
            }
        }
    }