I am trying to solve for X in BC=X where B is a symmetric positive semi-definite matrix and X is simply a vector with 1 in all of it’s entries. As far as I understand, this can be solved using C=B\X , but I was wondering that how can I take advantage of the positive semi-definite nature of B in order to make the solving faster or more numerically stable. I was using the below code

B = zeros(p, p)
for i in 1:p
for j in 1:i
B[i,j] = dot(R_iter_storage[i],R_iter_storage[j])
B[j,i] = B[i,j]
end
end
B = Symmetric(B)
X = ones(p)
C = B\X

I know that as I have declared B as a symmetric matrix, Julia would automatically use the optimized algorithm for symmetric matrices instead of the generalized solving algorithm. However, I have an additional constraint on B which is that it is a positive semi-definite matrix. Any ideas on what functions to use to solve the system such that the later property is also exploited ? Note - Here R_iter_storage is a vector whose elements are tensors of rank 4.

Please look into linearsolve. This package lets you solve system of linear equations very efficiently. Since your matrix is symm-positive definite, the method is recomended: LinearSolve.KrylovJL_CG()

Basically the command is :
using LinearAlgebra, LinearSolve, SparseArrays (if you have sparse vectors/matrices)
prob = LinearProblem(B, X);
sol = solve(prob,LinearSolve.KrylovJL_CG())

The matrices I would be working with are in the order of 6x6. Using the LinearSolve.KrylovJL_CG() seems to be causing numerical instability at the later stages of the DIIS solver for Coupled Cluster Equations I am trying to implement. The matrices I am trying to work with have small values of entries like

I would start with just cholesky(Hermitian(B)) \ X (you can use Symmetric instead of Hermitian, it makes no difference for real matrices). Using the Cholesky factorization will bestow the numerical and efficiency benefits you’re after.

Although if your matrix is singular (not strictly positive definite, or sufficiently close to singular that numerical issues arise) then cholesky will fail. In that situation, you will need to use pinv(B) * X or one of these iterative methods others have proposed.