symmetric positive definite). Updating $A^{p+1}_{ij} = A^p_{ij} U^{(p)}_{ij}$ orthogonalizes the block columns, since we have $${A^{p+1}_{ij}}^T A^{p+1}_{ij} = {U^{(p)}_{ij}}^T {A^p_{ij}}^T A^p_{ij} U^{(p)}_{ij} = {U^{(p)}_{ij}}^T G^{(p)}_{ij} U^{(p)}_{ij} = \Lambda^{p}_{ij},$$ where $\Lambda^{p}_{ij}$ is a diagonal matrix of the singular values of $G^{(p)}_{ij}$. Orthogonalizing all pairs of block columns until the entire matrix is orthogonal will give us the left singular vectors as the normalized columns and the singular values as the corresponding column norms. If the right singular vectors are needed, we can accumulate the action of the block rotation matrices on the identity matrix. For our batched implementation, we use highly optimized batched `syrk` and `gemm` routines from MAGMA to compute $G$ and to apply the block rotations, while the SVD is computed by our shared memory batched kernel. Since different matrices will converge in different numbers of sweeps, we keep track of the convergence of each operation $l$ by computing the norm $e_l$ of the off-diagonal entries of $G$ scaled by its diagonal entries. While this term is an inexact approximation of the off-diagonal terms of the full matrix in each sweep, it is still a good indication of convergence and will cost us at most an extra cheap sweep, 