def simulate_lu_decom(sim_locations,sample_locations,vmodel):
    c11 = fill_cova(sample_locations,None,vmodel)
    c21 = fill_cova(sim_locations,sample_locations,vmodel)
    c22 = fill_cova(sim_locations,None,vmodel)

    u11 = cholesky(c11)
    l11 = u11.T
    u11_inv = inv(u11)

    l21 = c21 @ u11_inv
    u12 = l21.T

    l22 = cholesky([email protected],lower=True)

    return u11_inv.T,l21,l22



    l11,u11 = lu(c11,permute_l= True)

    l11_inv = inv(l11)
    a21t = l11_inv @ c21.T
    a21 = a21t.T
    b12 = a21t

    l22,u22 = lu([email protected],permute_l= True)

    return a21,l11_inv,l22

def randomized_pca(A, n_components, n_oversamples=10, n_iter="auto",
                   flip_sign=True, random_state=0):
    """Compute the randomized PCA decomposition of a given matrix.

    This method differs from the scikit-learn implementation in that it supports
    and handles sparse matrices well.

    """
    if n_iter == "auto":
        # Checks if the number of iterations is explicitly specified
        # Adjust n_iter. 7 was found a good compromise for PCA. See sklearn #5299
        n_iter = 7 if n_components < .1 * min(A.shape) else 4

    n_samples, n_features = A.shape

    c = np.atleast_2d(A.mean(axis=0))

    if n_samples >= n_features:
        Q = random_state.normal(size=(n_features, n_components + n_oversamples))
        Q = safe_sparse_dot(A, Q) - safe_sparse_dot(c, Q)

        # Normalized power iterations
        for _ in range(n_iter):
            Q = safe_sparse_dot(A.T, Q) - safe_sparse_dot(c.T, Q.sum(axis=0)[None, :])
            Q, _ = lu(Q, permute_l=True)
            Q = safe_sparse_dot(A, Q) - safe_sparse_dot(c, Q)
            Q, _ = lu(Q, permute_l=True)

        Q, _ = qr(Q, mode="economic")

        QA = safe_sparse_dot(A.T, Q) - safe_sparse_dot(c.T, Q.sum(axis=0)[None, :])
        R, s, V = svd(QA.T, full_matrices=False)
        U = Q.dot(R)

    else:  # n_features > n_samples
        Q = random_state.normal(size=(n_samples, n_components + n_oversamples))
        Q = safe_sparse_dot(A.T, Q) - safe_sparse_dot(c.T, Q.sum(axis=0)[None, :])

        # Normalized power iterations
        for _ in range(n_iter):
            Q = safe_sparse_dot(A, Q) - safe_sparse_dot(c, Q)
            Q, _ = lu(Q, permute_l=True)
            Q = safe_sparse_dot(A.T, Q) - safe_sparse_dot(c.T, Q.sum(axis=0)[None, :])
            Q, _ = lu(Q, permute_l=True)

        Q, _ = qr(Q, mode="economic")

        QA = safe_sparse_dot(A, Q) - safe_sparse_dot(c, Q)
        U, s, R = svd(QA, full_matrices=False)
        V = R.dot(Q.T)

    if flip_sign:
        U, V = svd_flip(U, V)

    return U[:, :n_components], s[:n_components], V[:n_components, :]

def decompose( matrix ):
	# Returns the decomposition of a matrix A where
	#
	# Q.A.Q = P.L.U
	#
	# P.L.U is the factoring of Q.A.Q such that L is a lower triangular matrix with 1's
	# on the diagonal and U is an upper triangular matrix; P is the permutation (row-swapping
	# operations) required for this procedure. The permutation matrix Q is chosen such that 
	# the last element of U is its smallest diagnoal element. If A has a zero eigenvalue, 
	# then U's last element will be zero.
	
	dim = matrix.shape[ 0 ]

	# first decomposition
	( P, L, U ) = lu( matrix )
	
 	# detect the smallest element of U
	smallestIndex = findsmallestdiag( U )
	smallest = U[ smallestIndex, smallestIndex ]

	#show( matrix, "M" )
	#show( U, "U" )
	#print "Smallest element is %f at %d" % ( smallest, smallestIndex )

	# is the permutation Q not just the identity matrix?
	Q = identity( dim )
	if smallestIndex+1 != dim :
		# trick: exchange row 'smallestIndex' with row 'dim-1' of the identity matrix
		swaprow( Q, smallestIndex, dim-1 )

	return ( P, L, U, Q )

    def __get_Qd(self):
        """
        Sets the integration matrices QI and QE for the IMEX sweeper

        Returns:
            QI: St. Martin's trick, use LU decomposition
            QE: explicit Euler matrix, will also act on u0
        """
        # QI = np.zeros(np.shape(self.coll.Qmat))
        QE = np.zeros(np.shape(self.coll.Qmat))
        for m in range(self.coll.num_nodes + 1):
            # QI[m, 1:m+1] = self.coll.delta_m[0:m]
            QE[m, 0:m] = self.coll.delta_m[0:m]

        # This is for using LU decomposition
        # strip Qmat by initial value u0
        QT = self.coll.Qmat[1:,1:].T
        # do LU decomposition of QT
        [P,L,U] = LA.lu(QT,overwrite_a=True)
        # enrich QT by initial value u0
        Qd = np.zeros(np.shape(self.coll.Qmat))
        Qd[1:,1:] = U.T
        QI = Qd

        return QI, QE

def part_b(run_count, a, b):
    """
    Solve using LU decomposition
    """
    _, l, u = lu(a)
    for run in xrange(run_count):
        inv(u).dot(inv(l).dot(b))

    def __init__(self, num_nodes, tleft, tright):
        super(CollGaussRadau_Right_LU_Trick, self).__init__(num_nodes, tleft, tright)

        Q = self.Qmat

        p, l, u = lu(Q[1:, 1:].transpose())
        # print np.diag(l)
        self.QDmat = u.transpose()

def lu(A, b):
    #sol = []
    # Edit here to implement your code
    L,U=lp.lu(A,True)
    y=lp.solve(L,b)
    x=lp.solve(U,y)
   
    return list(x)

def find_basis(M):
    """Find the indices of the columns of M that form a basis or range(M)"""
    p,l,u = sla.lu(M)
    ind = [i for i in range(u.shape[0]) if u[i,i] != 0.0]
    if u[i,i] == 0:
        for j in range(i+1,u.shape[1]):
            if u[i,j] != 0.0: ind.append(j); break
    return ind

    def _traverse_grid_(self):
        """ Solve using linear systems of equations """
        P, L, U = linalg.lu(self.M1)

        for j in reversed(range(self.N)):
            x1 = linalg.solve(L, np.dot(self.M2, self.grid[1:self.M, j+1]))
            x2 = linalg.solve(U, x1)
            self.grid[1:self.M, j] = x2

 def _traverse_grid_(self):
     P,L,U=linalg.lu(self.coeffs)
     aux=np.zeros(self.M-1)
     
     for j in reversed(range(self.N)):
         aux[0]=np.dot(-self.a[1],self.grid[0,j])
         x1=linalg.solve(L,self.grid[1:self.M,j+1]+aux)
         x2=linalg.solve(U,x1)
         self.grid[1:self.M,j]=x2
         

 def __init__(self,Z,Az,r,apply='LU'):
     M=dgemm(Z,Az.T)
     if apply=='eig':
         self.setting_inverse_w_eigenvalues(M)
         super(CoarseLO,self).__init__(nargin=r,nargout=r,matvec=self.mult_eig,
                                         symmetric=True)
     elif apply =='LU':
         self.L,self.U=lu(M,permute_l=True,overwrite_a=True,check_finite=False)
         super(CoarseLO,self).__init__(nargin=r,nargout=r,matvec=self.mult,
                                         symmetric=True)

    def _traverse_grid_(self):
        """ Solve using linear systems of equations """
        P, L, U = linalg.lu(self.coeffs)
        aux = np.zeros(self.M-1)

        for j in reversed(range(self.N)):
            aux[0] = np.dot(-self.a[1], self.grid[0, j])
            x1 = linalg.solve(L, self.grid[1:self.M, j+1]+aux)
            x2 = linalg.solve(U, x1)
            self.grid[1:self.M, j] = x2

def Lu(A):
  """(L,U,P) = Lu(A)
  Compute pivoted LU decompostion of a matrix.

  RETURNS:
             L,U,P

      where A = PLU
  """
  (P,L,U) = lu(A)
  return ( L,U,P )

def lu(A, b):
    sol = []
    P, L, U = sci.lu(A)  
    #P is Permutation matrix 
    #L is Lower Triangular matrix
    #U is Upper Triangular matrix

    #A = PLU    
    #Ax = b  =>  PLUx = b
    sol = np.dot(sci.inv(U),np.dot(sci.inv(L),np.dot(sci.inv(P),b)))
    return list(sol)

def G2AH(G):
    n = len(G)
    A,B = npG2SVAR(G)
    P,L,U = linalg.lu(B)
    A = linalg.inv(L).tolist()
    B = B.tolist()
    A = listplace(A, 0.0, 0.0)
    for i in range(0,n): A[i][i] = 1
    B = listplace(B, 0.0, 'e')
    for i in range(0,n): B[i][i] = 'e'
    return A,B,P

def CompleteBase(V, B, eps=1e-4):
    tbase = append(V, B, axis=1)
    p, l, u = lu(tbase)
    echelon = zeros(u.shape[1], int)
    
    for row in u:
        tmp = nonzero(abs(row) > eps)[0]
        if tmp.size:
            echelon[tmp[0]] = 1
  
    return compress(echelon, tbase, axis=1)

def ldl_decomposition(matrix):
    p, L, U = scl.lu(matrix)
    size = matrix.shape[0]
    D = np.zeros((size, size), dtype=np.float64)
    divisor = np.ones(size, dtype=np.float64)
    for i in range(size):
        D[i][i] = U[i][i]
        divisor = U[i][i]
        U[i] /= divisor

    return L, D, U

def forward(A,b):
	p,l,u=s.lu(A)
	newB=b
	solutionList=[]
	for z in range(0,newB.shape[0]):
		mySum=0
		for j in range(0,z):
			mySum=mySum+ (l.item((z),j)*solutionList[j])
		xm=(newB.item(z)-mySum)/l.item(z,z)
		solutionList.append(xm)
	return mat(asarray(solutionList)).T

    def __refactorEtaFile(self):
        R = sps.identity(self.m)

        # multiply each eta-matrix
        for eta in self.etaFile:
            R = R * eta

        # decompose matrix
        R = R.todense()
        P, L, U = spl.lu(R)
        self.LU = [L, U]
        self.etaFile = []

def backwards(A,b):
	p,l,u=s.lu(A)
	newB=b
	solutionList=[]
	for z in range(0,newB.shape[0]):
		mySum=0
		for j in range(0,z):
			mySum=mySum+ (u.item(z,z-j)*solutionList[z-1-j])
		xm=(newB.item(newB.shape[0]-1 -z)-mySum)/u.item(newB.shape[0]-1 -z,newB.shape[0]-1 -z)
		solutionList.append(xm)
	solutionList=solutionList[::-1]
	return mat(asarray(solutionList)).T