What this equation is attempting to convey is that the matrix S ν kl µ ij is first computed from a slice of the global PAO overlap matrix S ν µ , and then transformed to the PNO basis through the PNO overlap matrices. In all our presented algorithms, every tensor operation that is not a dot product (i.e. two tensors sharing the same indices contracted to a scalar), can be performed through the level 3 BLAS matrix-multiplication subroutine GEMM.